Machine learning-based survival prediction using tumor lipid metabolism genes in osteosarcoma

Data collection and processing

We utilized the TCGAbiolinks package to retrieve gene expression and survival profiles of osteosarcoma samples from the TARGET database and explicitly selected the “TARGET-OS” project.¹³Additionally, we identified three independent validation datasets on osteosarcoma: GSE21257¹⁴GSE39058¹⁵and GSE16091¹⁶They were selected due to the availability of both gene expression and survival profiles. After this preprocessing step, the resulting number of samples was 84 for TARGET-OS, 53 for GSE21257, 37 for GSE39058, and 34 for GSE16091. When a gene had multiple values, we chose the median expression value to ensure data consistency.

Molecular subtype identification

Tumor lipid metabolism genes (TLMG) (from genecard)

In the Genecard database (https://www.genecards.org)[all] (Tumor) and [all] (Lipids and Metabolism) To identify genes related to lipid metabolism in tumors¹⁷From the Genecard database, we obtained a list of genes, each associated with a relevance score value. The relevance score is calculated by considering the importance of different resources that associate genes with diseases. To focus our study on genes that are strongly related to tumor lipid metabolism, we set a threshold and selected only genes with a relevance score above 10. This filtering criterion allowed us to identify a subset of genes, which we called tumor lipid metabolism genes (TLMGs). Subsequent analysis was limited to TLMGs with accessible gene expression data in the TARGET-OS, GSE21257, GSE39058, and GSE16091 datasets. After this initial screening, we found that 1,815 TLMGs had available gene expression values ​​across these datasets.

Molecular Subtypes of Osteosarcoma

To reveal distinct lipid metabolism patterns, we utilized unsupervised class detection by a consensus clustering algorithm facilitated by the “ConsensusClusterPlus” package.¹⁸This approach provides a solid foundation for determining the optimal number and configuration of clusters, increasing the accuracy of potential sample grouping. Notably, the sample sizes of the TARGET-OS, GSE21257, GSE39058, and GSE16091 datasets are relatively small, with the largest one containing only 84 samples. To increase the reliability of subtype identification, it was essential to combine samples from these four datasets into one unified dataset. However, directly combining these samples would pose a challenge called batch effect, which may distort the results. To mitigate this issue, we applied a binary transformation to the gene expression values ​​within each cohort, classifying them as 0 or 1 based on whether they were below or above the median expression level, respectively. This procedure, called binary transformation, has been proven to significantly minimize batch effect and is a strategy that has been widely adopted in research.¹⁹After this processing step, expression profiles from TARGET-OS, GSE21257, GSE39058, and GSE16091 datasets were aggregated into a unified meta-cohort. The clustering algorithm used the expression profiles of tumor lipid metabolism genes (TLMG) in the meta-cohort and employed the specified parameters (clusterAlg = “hc” for hierarchical clustering, distance = “pearson” for measuring correlation-based distance). To ensure the robustness of the classification, the clustering process was performed 500 times. This iteration allows us to validate the stability and consistency of the identified clusters and determine the optimal number of molecular subtypes required upon thorough evaluation of silhouette values, consensus matrices, delta areas, and tracking plots.

Differences in molecular subtypes

To clarify the variation in survival rates among different molecular subtypes, Kaplan-Meier survival curves supplemented with the log-rank test were used to evaluate the survival differences among these subtypes.Furthermore, we calculated enrichment scores for lipid metabolism pathways using the Gene Set Mutation Analysis (GSVA) ​​package.²⁰This method allowed for a detailed comparison of lipid metabolism scores among identified molecular subtypes, providing insight into the unique biological functions and pathways that define each subtype.In our previous study, 114 metabolic pathways were identified, of which 21 were dedicated to lipid metabolism, serving as a valuable resource for calculating lipid metabolism pathway enrichment scores.^{twenty one}Leveraging this information, we measured the activity levels of these 21 lipid metabolism pathways in the two molecular subtype samples from our study, which facilitated the comparison of lipid metabolism between the subtypes.

Using the pRRophetic package in R, we performed an extensive drug sensitivity analysis aimed at identifying differential sensitivity between the two subtypes of osteosarcoma. This analysis facilitated the identification of drugs that responded differently between these subtypes and helped select potential therapeutic agents tailored to specific osteosarcoma subtypes. Selection criteria were set rigorously, with emphasis on drugs that showed significant Benjamini and Hochberg adjusted sensitivity. p-values ​​less than 0.05 and significant group differences greater than 0.2.

Signature construction and verification

Feature Selection for Signatures

To address the issue of possible overfitting, we selected the TARGET-OS cohort from the four cohorts for feature/gene selection. First, we applied Fisher's exact test to analyze tumor lipid metabolism genes (TLMGs) between the two subtypes. TLMGs were selected because they are suitable for small binary data sets to identify genes with significant differences in expression. The results were ranked by p-value and adjusted using the Benjamini–Hochberg method to reduce false positives, focusing on genes with adjusted p-values ​​less than 0.05, indicating differential expression. These genes identified as differentially expressed genes (DEGs) were selected for further analysis. (2) Univariate Cox proportional hazards analysis was performed on the DEGs in the TARGET-OS cohort to evaluate their correlation with patient survival outcomes and identify genes with significant prognostic significance. (3) This method stepwise identified the most prognostic informative genes while emphasizing the accuracy and simplicity of the model. Through this process, we prioritized genes that conferred the greatest predictive value on the survival of osteosarcoma patients and identified and selected key genes for constructing a prognostic signature.

Construction of lipid metabolism-related signature (LMRS)

To develop a lipid metabolism-related signature (LMRS) with improved accuracy and robustness, we leveraged the strengths of four different machine learning algorithms to create 10 different algorithm configurations. The ensemble included four standalone algorithms, namely, random survival forest (RSF), CoxBoost, generalized boosted regression model (GBM), and survival support vector machine (SVM), as well as six combined approaches that integrated these techniques. The algorithm selection and combination was aimed at optimizing the predictive performance of the LMRS by leveraging the unique predictive capabilities and strengths of each method.

A comprehensive and sequential methodology was used to develop the signature. First, a combination of 10 different algorithms was applied to the prognostic genes selected in the step “2.3.1. Signature feature selection”. Following this first step, we thoroughly evaluated the performance of each algorithm across three different cohorts, specifically (1) TARGET-OS, (2) the combined cohort of GSE21257, GSE39058, and GSE16091, and (3) the meta-cohort. This approach ensured a thorough evaluation across different cohorts. Due to the relatively limited sample sizes of GSE21257, GSE39058, and GSE16091, it was necessary to combine them into one cohort. This combination was essential to accurately validate the performance of the lipid metabolism-related signature (LMRS), as performing 5-fold cross-validation on each dataset individually would have introduced errors due to lack of data. Notably, it was important to evaluate LMRS within the combined GSE cohort, as gene selection was finalized in the TARGET-OS cohort, with GSE21257, GSE39058, and GSE16091 providing completely independent datasets for testing LMRS. The final step in this process was to calculate the Harrell's concordance index (C-index) for each model across all validation datasets. The model exhibiting the highest C-index was then identified as optimal.

Evaluation of the LMRS model

The best model capable of generating predicted probabilities was used to define the LMRS score for each patient. Based on the median LMRS score, the LMRS scores of the patients in the different cohorts were calculated and then divided into high and low groups. This stratification validated the validity of the LMRS. To quantitatively evaluate the performance of the LMRS, the area under the receiver operating characteristic curve (AUC) values ​​were calculated. These values ​​measure the ability of the signature to effectively distinguish patients with different outcomes. Then, the log-rank test was performed to rigorously evaluate the difference in overall survival rates between the high and low LMRS groups. This statistical test allows comparing the survival distributions of the two groups and determining whether the survival times are significantly different. The analysis was complemented by the plot of survival curves, which visually represents the survival probability of both groups over time.

Comparison of LMRS with other publication signatures

To compare the performance of the lipid metabolism-related signatures (LMRS) with existing benchmarks, we conducted a thorough search of published signatures and ultimately obtained 12 distinct signatures. These signatures were developed using different algorithms, such as stepwise Cox regression and Lasso regression, reflecting a wide range of methodological approaches. Moreover, these signatures encompass a wide range of biological processes, including ferroptosis, stemness, fatty acid and lactate metabolism, aging, hexosamine biosynthesis, hypoxia, immune response, disulfidoptosis, drug sensitivity, and oxidative stress. This diversity highlights the multifaceted nature of the biological mechanisms they represent. For each of these signatures, we carefully calculated the Harrell's concordance index (C-index) across all cohorts to measure their predictive accuracy and efficacy in prognostic prediction. This comprehensive comparison highlights the relative performance of LMRS and demonstrates its potential advantages.

Relationship between LMRS score and molecular subtype

A T-test was used to compare the mean LMRS scores between the two molecular subtypes to identify significant differences. Box plots were then used to visually display the distribution of LMRS scores for each subtype. These plots display the median LMRS scores, variability, and outliers, providing a better understanding of how subtypes differ in terms of LMRS distribution.