Osteosarcoma (OS) is the most common primary malignant bone tumor, disproportionately affecting children and adolescents. The global incidence is approximately three cases per million people annually. Despite aggressive multimodal treatment combining neoadjuvant chemotherapy and surgery, the 5-year survival rate for patients who develop metastases remains dismal at only 10 to 20%. One of the core challenges in OS management is that patients with apparently identical clinical features, treatment regimens, and staging can have dramatically different survival outcomes. This molecular heterogeneity suggests that clinical variables alone are insufficient to guide prognosis and treatment intensification decisions.
The role of lipid metabolism: Lipids encompass a broad class of biomolecules including fatty acids, glycerides, phosphoglycerides, steroids, sphingolipids, and lipoproteins. In cancer biology, lipid metabolism has emerged as a critical driver of tumor progression, not merely a passive consequence of cellular growth. Lipid metabolism influences cancer cells through direct interactions with the tumor microenvironment and by activating oncogenic signaling pathways. Abnormal lipid accumulation suppresses dendritic cell function, impairing their ability to present tumor-associated antigens and activate anti-tumor T cells. In osteosarcoma specifically, diacylglycerols (a class of lipids) are overexpressed in metastatic compared to nonmetastatic cells, and inhibiting diacylglycerol synthesis reduces OS cell survival and motility. Lipid metabolism is measurably elevated in malignant and metastatic OS samples compared to less aggressive disease, making it a biologically plausible candidate for molecular subtyping.
The unmet need for molecular subtyping: High-throughput gene expression technologies now enable large-scale genomic profiling of tumors, but applying these data to define clinically actionable subtypes in osteosarcoma has been hampered by small dataset sizes and lack of multicenter validation. This 2024 study published in Scientific Reports addresses this gap directly: the authors use consensus clustering on lipid metabolism gene expression data pooled from four independent cohorts to define molecular subtypes, then build and rigorously validate a machine learning-based survival prediction signature derived from the subtype-discriminating genes.
The study aims to provide two deliverables for clinical translation: first, a biologically coherent molecular subtype classification of OS based on lipid metabolic programs, and second, a compact 12-gene prognostic signature that can stratify patients into high-risk and low-risk groups with strong predictive accuracy across independent validation cohorts.
The study assembled four publicly available osteosarcoma datasets: TARGET-OS (84 samples from the Therapeutically Applicable Research to Generate Effective Treatments database), GSE21257 (53 samples), GSE39058 (37 samples), and GSE16091 (34 samples). All four datasets included both gene expression profiles and survival outcomes, allowing consistent analysis. The combined Meta-Cohort contained 208 osteosarcoma samples, making it among the larger consolidated OS genomic datasets used for this type of analysis, given the rarity of the disease. Only samples with complete expression and survival data were retained after preprocessing.
Identifying Tumor Lipid Metabolism Genes (TLMGs): The authors queried the GeneCards database using the search term "[all] (tumor) AND [all] (lipid AND metabolism)" to compile a list of genes linked to lipid metabolism in the cancer context. Each gene in GeneCards carries a relevance score reflecting the strength of its association with disease-related resources. By setting a threshold of relevance score greater than 10, they identified 2,759 high-confidence TLMGs. After filtering to retain only those with available expression data across all four OS cohorts, 1,815 TLMGs were carried forward for downstream analysis. This gene set anchors the entire study in lipid biology rather than generic transcriptomic variation.
Handling batch effects with binary transformation: Directly merging gene expression data from four studies performed at different institutions and on different platforms introduces systematic batch effects that can generate spurious clusters. To address this, the authors applied binary transformation to each cohort's expression data, coding each gene as 0 (below cohort median) or 1 (above cohort median). This standardization removes absolute expression level differences driven by platform or laboratory variation, retaining only relative expression patterns. Binary transformation is a validated strategy for cross-cohort integration in studies with small, heterogeneous datasets. The four datasets were then merged into a single Meta-Cohort expression matrix of 1,815 TLMGs across 208 samples for subtype discovery.
Consensus clustering setup: Unsupervised molecular subtype discovery was performed using the ConsensusClusterPlus package in R, applying hierarchical clustering (clusterAlg = "hc") with Pearson correlation-based distance. The algorithm was run 500 times to ensure cluster stability. Optimal cluster number was evaluated using silhouette width plots, consensus cumulative distribution function (CDF) curves, delta area plots, and tracking plots. These multiple evaluation criteria were used in concert because different metrics can suggest different optimal cluster numbers in small datasets.
Consensus clustering of the 208-sample Meta-Cohort using 1,815 TLMG expression profiles identified two optimal molecular subtypes, labeled C1 and C2. The tracking plot and silhouette width both confirmed two as the most stable partition, even though the consensus CDF and delta area plots pointed toward four as a candidate optimal number. The investigators deliberately chose two subtypes based on two pragmatic considerations: the tracking plot showed that even when four clusters were explored, the vast majority of samples still segregated into two dominant groups, and a smaller number of subtypes simplifies downstream machine learning model construction and clinical application. In the Meta-Cohort, 115 samples were assigned to C1 and 93 to C2.
Survival differences: Kaplan-Meier survival analysis with log-rank testing showed that patients in C2 had significantly longer overall survival than those in C1 (p-value not specified beyond significance threshold). This survival divergence validates that the two lipid metabolism-based subtypes capture biologically meaningful prognostic differences rather than arbitrary statistical partitions.
Pathway-level differences between C1 and C2: Gene Set Variation Analysis (GSVA) was applied to calculate enrichment scores for 21 lipid metabolism pathways across both subtypes. Of the 21 pathways, 14 showed statistically significant differences between C1 and C2 by t-test. The C1 subtype is defined by elevated cholesterol metabolism, fatty acid biosynthesis and elongation, and ketone metabolism. These are energy-generating and membrane-building lipid programs typically associated with rapidly proliferating cells. The C2 subtype, in contrast, shows predominance of steroid hormone biosynthesis (including aldosterone, cortisol, estradiol, and testosterone), arachidonic acid metabolism, and glycerolipid and linoleic acid metabolism. These programs are more characteristic of differentiated cellular states and immune-modulatory lipid signaling, suggesting that C2 tumors occupy a different metabolic niche with distinct therapeutic vulnerabilities.
Drug sensitivity differences: Using the pRRophetic package, which applies ridge regression to estimate IC50 values from GDSC drug sensitivity data, the authors quantified differential drug sensitivity between the two subtypes. The analysis found that 48 drugs showed significantly higher sensitivity in the C1 subtype (lower IC50, adjusted p-value less than 0.05, group difference greater than 0.2), while only 5 drugs were more effective against C2. This asymmetry suggests that C1, despite having worse survival, may be more broadly targetable with currently available compounds, potentially informing future subtype-stratified clinical trials.
With 1,815 TLMGs as starting candidates, building a reliable prognostic signature required systematic dimensionality reduction to avoid overfitting. The authors restricted the entire feature selection pipeline to the TARGET-OS cohort (84 samples) exclusively, reserving the other three cohorts as independent validation sets. This design choice is critical: using separate data for feature selection and validation prevents the optimistic bias that arises when validation sets overlap with training data. The feature selection proceeded through three sequential stages, each progressively narrowing the gene list.
Stage 1 - Differential expression between subtypes: Fisher's exact test was applied to all 1,815 TLMGs in TARGET-OS to identify genes with significantly different binary expression profiles between C1 and C2 subtypes. Results were adjusted using the Benjamini-Hochberg method to control the false discovery rate. Genes with adjusted p-value below 0.05 were retained as differentially expressed genes (DEGs). This step yielded 698 genes that were meaningfully different between the two lipid metabolism subtypes within the discovery cohort.
Stage 2 - Univariate Cox analysis for prognostic relevance: The 698 DEGs were subjected to univariate Cox proportional hazards analysis in TARGET-OS to identify those with individual association with overall survival. This reduced the list to 35 genes with prognostic significance. Of these, 17 had negative Cox coefficients (protective, higher expression associated with better survival) and 18 had positive coefficients (risky, higher expression associated with worse survival). This directional classification is important for understanding the functional role of each gene in the signature.
Stage 3 - Stepwise AIC for model parsimony: The 35 prognostically significant genes were further refined using stepwise Akaike Information Criterion (stepAIC), an iterative variable selection procedure that adds or removes genes based on their impact on the AIC score. StepAIC balances model complexity against explanatory power, penalizing models with excess variables. The result was a final set of 12 genes that provided the optimal balance of prognostic information and model simplicity. These 12 genes served as the input features for machine learning signature construction.
The Lipid Metabolism-Related Signature (LMRS) was constructed using a methodologically rigorous ensemble approach that tested ten different machine learning algorithm combinations to identify the most predictively robust model. The four base algorithms used were: Random Survival Forest (RSF), CoxBoost (a gradient boosting method applied to Cox regression), Generalized Boosted Regression Models (GBM), and Survival Support Vector Machine (survival-SVM). These four algorithms were used both individually and in six pairwise combinations, yielding ten total configurations. Each combination was evaluated within a five-fold cross-validation framework applied to three cohort settings: TARGET-OS (84 samples), the combined independent validation cohort of GSE21257 plus GSE39058 plus GSE16091 (124 samples), and the full Meta-Cohort (208 samples).
Why combine cohorts for validation: The three independent validation cohorts (53, 37, and 34 samples respectively) were too small to support reliable five-fold cross-validation individually, as splitting 37 or 34 samples into five folds would leave fewer than 8 patients per fold, causing unstable training and testing estimates. Combining them into a single 124-sample independent validation cohort allowed robust cross-validation while preserving their independence from the TARGET-OS discovery cohort, since gene selection had already been finalized before these data were touched.
Model selection criterion: All ten algorithm combinations were assessed using Harrell's Concordance Index (C-index), which measures the proportion of patient pairs for which the model correctly ranks relative survival times. A C-index of 1.0 indicates perfect discrimination; 0.5 equals random prediction. The selection criterion was the highest average C-index in the independent testing cohorts (not the discovery cohort), specifically to prioritize generalizability over in-sample fit and to counteract overfitting.
Winning model: The combination of randomForestSRC and CoxBoost achieved the highest average C-index of 0.713 across the independent validation cohorts and was selected as the optimal LMRS model. For comparison, the CoxBoost plus GBM combination achieved a C-index of 0.874 within TARGET-OS, but this dropped to 0.674 in the independent validation cohorts, illustrating exactly the overfitting problem the selection criterion was designed to avoid. The randomForestSRC and CoxBoost model sacrificed some in-sample performance for substantially better generalization.
With the optimal randomForestSRC plus CoxBoost model established, each patient across all cohorts received a continuous LMRS risk score. Patients were dichotomized into high-score and low-score groups based on the median LMRS value in each respective cohort. Kaplan-Meier analysis consistently showed that patients in the high-score group had significantly worse overall survival than those in the low-score group, with statistical significance (all p-values less than 0.05) across the Meta-Cohort, TARGET-OS, and the combined GSE validation cohort. This consistent separation across three independent datasets demonstrates that the signature captures prognostic signal that transfers beyond its discovery context.
Time-dependent ROC analysis: The predictive accuracy of LMRS was quantified using area under the receiver operating characteristic curve (AUC) calculated at 1, 3, and 5 years. In the Meta-Cohort, AUC values were 0.745, 0.743, and 0.730 at these time points respectively. In TARGET-OS, the AUCs were 0.765, 0.788, and 0.785. In the combined GSE validation cohort, they were 0.741, 0.708, and 0.692. The AUC values show a consistent pattern: strong discrimination at all three survival horizons across all cohorts, with relatively small drops from internal to external validation, indicating the signature generalizes well despite the small sample sizes.
Benchmark against 12 published signatures: A rigorous head-to-head comparison was conducted against 12 previously published osteosarcoma prognostic signatures covering diverse biological processes including ferroptosis, cellular stemness, fatty acid and lactate metabolism, aging, hexosamine biosynthesis, hypoxia, immune response, disulfidptosis, drug sensitivity, and oxidative stress. For fair comparison, each published signature was evaluated using the same ten machine learning algorithm combinations and the same three cohort sets. LMRS achieved the highest C-index values across all three settings: 0.748 in the Meta-Cohort, 0.874 in TARGET-OS, and 0.713 in the combined GSE cohort. In comparison, the best-performing published signature reached maximum C-indices of 0.664 (Meta-Cohort), 0.815 (TARGET-OS), and 0.700 (combined GSE). LMRS outperformed all 12 competitors in the Meta-Cohort and combined GSE cohort, and came second to one in TARGET-OS.
Gene-level risk contributions: Univariate Cox regression of the 12 signature genes in the Meta-Cohort identified directional risk contributions. MUC1 (mucin 1), PPP2R1B (protein phosphatase 2 regulatory subunit A beta), and ME1 (malic enzyme 1) were associated with worse survival when expressed at higher levels. ABCD3 (ATP-binding cassette transporter D3) and CD4 (cluster of differentiation 4, a T-cell marker) were associated with better survival at higher expression, potentially reflecting favorable immune cell infiltration patterns in the tumor microenvironment.
An important internal consistency check for any molecular classification-plus-signature workflow is whether the derived signature scores actually align with the original subtype labels. If the LMRS scores are independent of the C1/C2 subtype assignment, that would suggest the signature captures different biology than the subtype clustering. By t-test comparison across the Meta-Cohort, the C1 subtype had significantly higher LMRS scores than C2. This is directionally coherent: C1 was identified as having worse survival in the Kaplan-Meier analysis, and high LMRS scores also correspond to worse survival in subsequent validation. The two analytical frameworks therefore converge on the same biological signal, providing mutual validation.
Biological interpretation: The C1 subtype's elevated LMRS scores reflect its association with cholesterol, fatty acid biosynthesis, and ketone metabolism programs. Within the 12-gene LMRS, elevated MUC1 expression is consistent with aggressive tumor biology, as MUC1 is a transmembrane glycoprotein that promotes oncogenic signaling and suppresses immune recognition across multiple cancer types. ME1 catalyzes the conversion of malate to pyruvate with NADPH production, linking lipid biosynthesis to central carbon metabolism. PPP2R1B is a regulatory subunit of protein phosphatase 2A (PP2A), and its loss-of-function role as a tumor suppressor in various cancers makes elevated expression counterintuitive, though context-dependent PP2A activity can vary by subtype and phosphoproteome state.
Conversely, the protective roles of ABCD3 and CD4 in LMRS are biologically interpretable. ABCD3 is a peroxisomal membrane transporter involved in fatty acid oxidation and is part of a metabolic axis that can antagonize the lipid biosynthesis programs enriched in C1. CD4's protective association may reflect CD4-positive helper T-cell infiltration in the tumor microenvironment, which is a recognized favorable prognostic feature in several solid tumors including osteosarcoma.
The alignment between LMRS scores and molecular subtypes demonstrates that the signature is not merely a statistical abstraction but reflects the same underlying lipid metabolic programs that define the biological subtypes, strengthening confidence in its mechanistic basis and clinical translatability.
Small and geographically restricted cohorts: Osteosarcoma is a rare disease with an incidence of approximately 3 cases per million people annually. This rarity fundamentally limits the scale of publicly available genomic datasets. The largest single cohort used in this study (TARGET-OS) contains only 84 samples. Even the combined Meta-Cohort of 208 samples is modest by contemporary machine learning standards. While the binary transformation and multicenter design partially mitigate batch effects and promote generalizability, the small sample sizes increase the variance of performance estimates and reduce statistical power for detecting subgroup effects. The three independent validation cohorts (37, 53, and 34 samples) had to be pooled to support five-fold cross-validation, which means even the "independent" validation partially reduces geographic independence.
Lack of experimental validation: The 12 LMRS genes were selected purely through computational analysis of gene expression data. No functional experiments were performed to confirm that any of these genes causally drive the lipid metabolic differences between C1 and C2, or that perturbing their expression changes osteosarcoma cell behavior. The paper explicitly acknowledges that experimental validation (cell line knockdown or overexpression experiments, animal models, immunohistochemistry in tumor tissue) is required to support the computational findings before any clinical application. Without this, the signature remains correlative rather than mechanistically confirmed.
Geographic and racial generalizability: The TARGET-OS cohort draws primarily from Norwegian and American patients, while the GEO datasets also represent predominantly Western patient populations. Osteosarcoma biology may differ across geographic and ethnic groups due to differences in genetic background, environmental exposures, and treatment protocols. The authors specifically flag that application of the LMRS to Chinese patients, for example, is uncertain, and plan to address this through future multicenter collaborations that include Chinese patient cohorts. This is a critical validation step before any pan-ethnic clinical deployment.
Future directions: The authors propose several next steps. First, multicenter validation in Chinese and other Asian patient populations to test geographic generalizability. Second, wet-lab experiments to functionally characterize the 12 LMRS genes and establish causal roles in lipid metabolism and OS progression. Third, integration of the LMRS with clinical variables (stage, metastasis status, treatment response) into a multivariate nomogram that could be tested in prospective settings. Fourth, exploration of whether the two lipid metabolism subtypes respond differently to emerging lipid pathway inhibitors, given the drug sensitivity differences (48 C1-sensitive vs. 5 C2-sensitive compounds) identified computationally.