Prognostic value of a glycolysis and cholesterol synthesis related gene signature in osteosarcoma: implications for immune microenvironment and personalized treatment strategies

Data collection and processing

The osteosarcoma cohort (TARGET-OS) was acquired from The Cancer Genome Atlas (TCGA) database (https://portal.gdc.cancer.gov/), encompassing 88 tissue samples. The obtained data were in the form of Fragments Per Kilobase of exon model per Million mapped fragments (FPKM) values, which underwent logarithmic transformation for subsequent analysis.

To further enrich our investigation, we gathered a collection of 72 genes associated with “REACTOME_GLYCOLYSIS” and 27 genes associated with “REACTOME_ CHOLESTEROL _BIOSYNTHESIS” from the Gene Set Enrichment Analysis (GSEA) website available at http://www.gsea-msigdb.org. Upon identification of one undetected gene within the TCGA cohort, our final gene set consisted of 98 genes, forming the basis of the glycolysis and cholesterol biosynthesis-related gene set for subsequent analysis.

Unsupervised clustering analysis of TARGET-OS samples

To explore the heterogeneity within the 88 TARGET-OS samples, we performed an extensive unsupervised clustering analysis using the Consensus Clustering method. This approach significantly enhances the precision and reliability of deciphering complex patterns within the dataset [23, 24]. The assessment of the optimal number of subgroups was based on evaluating the Cumulative Distribution Function (CDF) values and carefully inspecting the clustering heatmap.

Following the identification of distinct subgroups, we conducted a comprehensive differential gene expression analysis by considering all genes associated with osteosarcoma. Additionally, we compared the overall survival rates between these subgroups. To specifically investigate the differential expression patterns of the 98 GCSRG genes within the identified subgroups, we employed the “limma” R package to perform a differential analysis.

Analysis of mutation frequency and tumor mutation burden in different subgroups

To gain a comprehensive understanding of the genetic characteristics underlying osteosarcoma and identify potential therapeutic targets, we conducted an in-depth analysis of mutation frequencies within the distinct GCSRG subgroups. By leveraging the rich genomic data available from the TCGA database, our objective was to uncover the divergent mutational profiles between these subgroups and their potential implications for clinical outcomes.

Furthermore, we conducted an extensive examination of potential variations in tumor mutation burden (TMB) between the two subgroups using the “maftools” R package. This comprehensive analysis allowed us to evaluate the significance of TMB in relation to osteosarcoma, shedding light on its potential as a prognostic biomarker and its influence on the immune microenvironment and response to treatment.

Tumor Immune Dysfunction and Exclusion (TIDE) score analysis

We scrutinized the immune therapy response between the high- and low-risk groups by analyzing the Tumor Immune Dysfunction and Exclusion (TIDE) (http://tide.dfci.harvard.edu/) scores, which could provide indications of immunotherapy responsiveness in osteosarcoma patients.

Assessment of immune microenvironment in different subtypes

In order to investigate the association between the two subtypes of TARGET-OS samples and immune cell subpopulations, we utilized the CIBERSORT algorithm. This approach enabled us to comprehensively evaluate the infiltration of 22 distinct immune cell types within the TARGET-OS samples. Furthermore, we conducted an in-depth analysis of the expression patterns of human leukocyte antigen (HLA)-related genes and immune checkpoint-related genes across the different subtypes, providing valuable insights into their immunological characteristics.

Regression model construction for GCSRG-related genes

Our regression model for the GCSRG-related genes was meticulously developed through a series of comprehensive steps. Initially, a univariate Cox regression analysis was performed to identify genes significantly associated with overall survival, utilizing a significance threshold of p-value <0.05. Subsequently, we refined the selection of survival-related genes by employing LASSO regression and multivariate Cox regression analyses. These rigorous techniques enabled us to construct a robust prognostic model that accurately estimates the prognosis of osteosarcoma patients based on the expression levels of genes [25]. This meticulous approach ensured the inclusion of only the most crucial genes, enhancing the reliability and clinical utility of our prognostic model.

GSEA enrichment analysis

To elucidate the distinct molecular processes that differentiate the high- and low-risk groups, we conducted Gene Set Enrichment Analysis (GSEA). This approach enabled us to identify enriched signaling pathways within each group, providing valuable insights into the underlying biological mechanisms. By employing GSEA, we sought to unravel the intricate molecular landscape and shed light on the key pathways driving the divergent risk profiles observed in our study.

Statistical processing

All statistical analyses were conducted using R software, version 4.1.0, developed by the R Foundation for Statistical Computing, Vienna, Austria. The Student’s t-test or Wilcoxon test was employed to compare the two groups. The Kruskal–Wallis test was used to analyze comparisons between several groups. The log-rank test and Kaplan–Meier plots, which display survival curves, were compared. If the p-value was less than 0.05, statistical findings were declared significant.

Optimal subgroup identification

Based on the CDF values and clustering heatmap, it was determined that the optimal number of subgroups for our analysis was two. This partitioning maximized the dissimilarities between the groups while minimizing variations within each group (Figure 1A and B).

Figure 1:

Identification of optimal subgroups based on gene expression profiles. (A) Clustering heatmap showing the gene expression pattern of the identified subgroups. (B) Cumulative distribution function (CDF) plot for determining the optimal number of subgroups. (C) Heatmap illustrating the differential gene expression between the two subgroups. (D) Kaplan–Meier survival curve comparing overall survival rates of the GCSRG low and high subgroups.

Through differential gene expression analysis, we identified substantial alterations in gene expression patterns between the two groups (Figure 1C). Additionally, when comparing the overall survival rates, we observed a significantly higher overall survival in the GCSRG low group (Figure 1D).

Mutation frequency analysis in different subgroups and TIDE score analysis

The analysis unveiled a notable disparity in mutation frequency between the GCSRG low subgroup (56 %) and the GCSRG high subgroup (80 %). TP53 and MUC16 emerged as the two most frequently altered genes across both subgroups, underscoring their potential significance (Figure 2A and B).

Figure 2:

Mutation frequency analysis in the GCSRG low and high subgroups. (A) Mutation frequency plot for the GCSRG low subgroup. (B) Mutation frequency plot for the GCSRG high subgroup. (C) Box plot comparing the tumor mutation burden between the two subgroups. (D) Box plot comparing the TIDE scores between the two subgroups.

Contrary to our expectations, no statistically significant distinction in tumor mutation burden was observed when comparing the two groups (Figure 2C). However, intriguingly, the high-risk group exhibited lower TIDE scores, indicative of a higher likelihood of favorable response to immunotherapy (Figure 2D).

Immune microenvironment in different subtypes

The bar plot in Figure 3A effectively visualized the abundance and distribution of the 22 distinct immune cell types within the TARGET-OS samples. Subsequent correlation heatmap analysis unveiled predominantly negative correlations among the various immune cell populations (Figure 3B). Notably, employing the ESTIMATE algorithm, we observed elevated stromal scores, immune scores, and microenvironment scores in the GCSRG low subtype, indicating a more dynamic and active microenvironment in this subgroup (Figure 3C–F). Furthermore, our assessment of HLA-related genes revealed higher expression levels in the GCSRG low subtype compared to the GCSRG high subtype (Figure 3G).

Figure 3:

Immune microenvironment analysis of the different subtypes. (A) Bar plot illustrating the infiltration of 22 immune cell types in the TARGET-OS samples. (B) Correlation heatmap of immune cell infiltration. (C–F) Box plots showing stromal scores, immune scores, and microenvironment scores determined using the ESTIMATE algorithm. (G) Box plot displaying the expression levels of HLA-related genes in the GCSRG low and high subtypes. (H) Box plot showing the expression levels of immune checkpoint-related genes in the two subtypes. *Indicates a p-value less than 0.05, **denotes a p-value less than 0.01, and ***signifies a p-value less than 0.001.

Consistent with our previous findings, the differential analysis of immune checkpoint-related genes demonstrated heightened expression levels in the GCSRG low subtype (Figure 3H). This observation further supports the notion that the GCSRG low subtype exhibits a heightened susceptibility to immune cell infiltration. Consequently, patients harboring this particular subtype may potentially benefit from enhanced effectiveness of immunotherapy interventions.

Regression model construction for GCSRG-related genes

Through the utilization of univariate Cox regression analysis, we successfully identified eight genes that displayed a significant association with overall survival (Figure 4A). To refine our selection of survival-related genes, we employed two filtering approaches. The first involved LASSO regression, where eight genes were retained based on the optimal lambda value, thus establishing a more precise gene subset (Figure 4B and C). Subsequently, employing multivariate Cox regression analysis, we established a prognostic model comprising eight genes: GNPDA2, NUP98, PFKFB2, PFKFB3, PKLR, PRKACB, SQLE, and TM7SF2.

Figure 4:

Regression model construction for GCSRG-related genes. (A) Forest plot of univariate Cox regression analysis identifying overall survival-related genes. (B–C) LASSO regression plot for selecting survival-related genes. (D) Kaplan–Meier survival curve comparing overall survival rates between high- and low-risk groups based on the risk scores. (E) ROC curve for evaluating the sensitivity and accuracy of the risk scores. (F–G) Plots illustrating the distribution of risk scores and survival status. (H) Heatmap of the expression levels of the eight genes in the high- and low-risk groups.

Risk scores, derived from the expression levels of the aforementioned eight genes and their respective regression coefficients, were calculated to stratify patients into high- and low-risk groups. Notably, these risk groups exhibited substantial disparities in overall survival outcomes (Figure 4D). Furthermore, the ROC curve analysis demonstrated that the risk score displayed remarkable sensitivity and accuracy in distinguishing between high- and low-risk patients (Figure 4E).

By utilizing the risk scores, we implemented a ranking system for patients and observed a consistent pattern where higher risk scores were associated with decreased survival rates (Figure 4F and G).

This aligns with the expected characteristics of a risk model. Among the eight genes included in the model, TM7SF2, SQLE, PKLR, and PFKFB3 exhibited high expression levels in the high-risk group, while the remaining four genes showed lower expression levels in the high-risk group (Figure 4H).

Prognostic model construction

In our study, univariate Cox regression analysis revealed that metastasis and risk score achieved statistical significance with a p-value below 0.05 (Figure 5A). Subsequently, in the multivariate Cox regression analysis, both metastasis and risk score emerged as independent prognostic indicators, underscoring their importance in predicting patient outcomes (Figure 5B). Notably, when assessing the performance of the four components through ROC curve analysis, the risk score exhibited the highest AUC value of 0.85 (Figure 5C).

Figure 5:

Prognostic model construction and TIDE score analysis. (A) Forest plot of univariate Cox regression analysis identifying factors with p-value <0.05. (B) Forest plot of multivariate Cox regression analysis determining independent prognostic factors. (C) ROC curve analysis of the four factors. (D) Construction of the prognostic prediction model based on the four factors. (E) Calibration curve of the prognostic prediction model. (F) Harrell’s C-index value calculation. (G) Box plot showing TIDE scores for the high- and low-risk groups, *** denotes a p-value less than 0.001.

Using the aforementioned four parameters, we constructed a comprehensive prognostic prediction model (Figure 5D). The calibration curve analysis further demonstrated excellent discrimination and calibration of the model (Figure 5E). Remarkably, when computing Harrell’s C-index, the risk score emerged as the parameter with the highest value, underscoring its superior predictive capability (Figure 5F).

GSEA enrichment analysis

The GSEA analysis revealed distinct enrichment patterns in signaling pathways between the low- and high-risk groups. Specifically, the low-risk group exhibited pronounced enrichment in signaling pathways such as the Jak–STAT signaling pathway, cytokine-cytokine receptor interaction, and neuroactive ligand–receptor interaction (Figure 6A). In contrast, the high-risk group demonstrated a prominent enrichment of pathways associated with Huntington’s disease, olfactory transduction, and other relevant pathways (Figure 6B).

Figure 6:

GSEA enrichment analysis results. (A) Enrichment plot for signaling pathways enriched in the low-risk group. (B) Enrichment plot for signaling pathways enriched in the high-risk group.