Construction of an oxidative phosphorylation-related gene signature for predicting prognosis and identifying immune infiltration in osteosarcoma

Training and validation data sets of osteosarcoma patients

The flowchart of this article is shown in Supplementary Figure 14. The current investigation merged the clinical and transcriptional data of 88 individuals with osteosarcoma retrieved from the UCSC Xena database. A total of 84 individuals diagnosed with osteosarcoma, along with their clinical information, were included, after excluding samples with unreported or zero survival time. The dataset of GSE21257 contained 53 cases, and the gene expression matrix obtained from the combination of clinical data and transcriptional data was used as the validation set. The clinical data of the two groups are shown in Table 1.

Prognosis-related oxidative phosphorylation genes and risk-scoring model of osteosarcoma

The GSE28425 dataset was used to obtain 3662 differential genes. The intersection with 342 oxidative phosphorylation-related genes was taken to obtain 68 oxidative phosphorylation-related genes (Figure 1A–1C). A univariate Cox regression analysis was conducted on the training set of individuals to identify five oxidative phosphorylation genes that were associated with prognosis: ATP6V0D1, LHPP, COX6A2, MTHFD2, and NDUFB9 (Figure 2B). Subsequently, the risk-scoring model was developed based on LASSO regression analysis (Figure 2C, 2D). Therefore, the present research computed the risk score of all patients. The correlation coefficients between the five oxidative phosphorylation-related differential expression prognostic genes and risk scores were also examined (Supplementary Figures 1, 2). Analysis of the TCGA training group data demonstrated that individuals with osteosarcoma classified as high-risk exhibited a remarkably poorer prognosis in comparison to those categorized as low-risk (P = 6.602e-03, Figure 3A). Scatter plots and risk curves showed survival and risk scores for all osteosarcoma samples. Individuals with osteosarcoma in the low-risk group demonstrated lower mortality rates and risk coefficients in contrast with the other group (Figure 3B, 3C). Moreover, the heat map highlighted the expression of five oxidative phosphorylation genes related to prognosis in 84 osteosarcoma patients in both groups (Figure 3D). A time-dependent ROC curve was utilized to compute the risk score, yielding respective AUC values of 0.749, 0.773, and 0.744 to predict the OS of individuals suffering from osteosarcoma at 1, 3, and 5 years (Figure 3E). PCA dimensionality reduction analysis showed that the risk-scoring model can better classify osteosarcoma patients into two groups, with considerable variations across the two groups (Figure 3F). The validation set in GSE21257 also demonstrated the effectiveness of the risk model to predict the survival of individuals suffering from osteosarcoma, with low-risk individuals exhibiting a longer OS period in contrast with low-risk individuals (P = 1.357e-02, Figure 4A). Additionally, scatter plots further illustrated that individuals at higher risk demonstrated a more unfavorable prognosis (Figure 4B). Risk curve plots, heat maps, and PCA dimensionality reduction analysis charts showed similar results to the training set (Figure 4C, 4D, 4F). Furthermore, the respective AUC values predicted by the validation group time ROC curve for the OS of individuals suffering from osteosarcoma for 1-, 3-, and 5-year periods were 0.735, 0.760, and 0.731 (Figure 4E). The risk-scoring model demonstrated excellent predictive performance, as evidenced by the AUC values exceeding 0.7 for both validation and training groups.

Figure 1. Acquisition of oxidative phosphorylation-related genes. (A, B) GSE28425 data set: Analysis of differences between osteosarcoma and adjacent cancer (osteoblasts) (differential gene volcano map and differential gene heat map). (C) Obtaining a Wayne map of differentially expressed genes related to oxidative phosphorylation in osteosarcoma.

Figure 2. Construction of a risk-scoring model for patients with osteosarcoma based on oxidative phosphorylation genes. (A) Forest plots based on univariate Cox regression analysis showed that there were 5 oxidative phosphorylation genes with prognostic significance (p < 0.01). (B, C) LASSO regression to construct a risk model.

Figure 3. TCGA osteosarcoma data were used as training groups to validate the effectiveness of the risk score model in predicting survival. (A) Kaplan-Meier survival analysis; (B) Scatter chart: green represents survival during follow-up, red represents death during follow-up, abscissa represents risk score, and ordinate represents survival time. (C) Risk curve. (D) Heat map of gene expression in 84 patients with osteosarcoma. (E) Time ROC graph. (F) PCA dimensionality reduction analysis chart.

Figure 4. GSE21257 data was used as a validation group to verify the effectiveness of the risk scoring model in predicting survival. (A) Kaplan-Meier survival analysis; (B) Scatter chart: green represents survival during follow-up, red represents death during follow-up, abscissa represents risk score, and ordinate represents survival time. (C) Risk curve. (D) Heat map of gene expression in 53 patients with osteosarcoma. (E) Time ROC graph. (F) PCA dimensionality reduction analysis chart.

Subgroup analysis

The prognosis-predictive efficacy of the risk-scoring model in various clinical feature subgroups was elucidated through the survival analysis of clinical features, with the following findings: Male (p = 0.036), female (p = 0.076), non-metastatic (p = 0.052), age 18 (p = 0.068), age 18 (p = 0.023), and metastatic (p = 0.050) (Supplementary Figure 3A–3F). The validation set application also involved conducting survival analysis based on clinical features, which elucidated the capability of the risk-scoring model to predict outcomes in different clinical feature subgroups (Supplementary Figure 4A–4F). Although P-values were not less than 0.05 for many subgroups, it was attributed to the limited sample size of patients and the apparent poorer survival observed among individuals in the high-risk group.

Tumor microenvironment across high- and low-risk groups

In order to explore the underlying reasons behind the significant impact of oxidative phosphorylation-related genes on the prognosis of individuals with osteosarcoma, the current research analyzed the results of the infiltration of 22 immune cells across the two risk groups. The results demonstrated the proportion of immune cell composition in 84 patients, highlighting that macrophages accounted for a large proportion of osteosarcoma samples. Macrophages have been found to be substantially involved in the onset and progression of osteosarcoma (Figure 5A). Moreover, this study identified five distinct types of immune cells that exhibited variations across the two groups. In the high-risk group, there was a higher infiltration level of CD4 naive T cells while a lower infiltration level of CD8+T cells, follicular helper T cells, activated CD4 memory T cells, and γ δ T cells was noted (Figure 5B). Subsequently, the link between immune cells was compared, and it was found that M2 macrophages were inversely proportional to most immune cells in osteosarcoma. The correlation coefficient with M0 macrophages was -0.58, showing a high negative correlation, while the correlation coefficient with M1 macrophages was 0.08, showing no correlation. The correlation coefficient between M1 macrophages and M0 macrophages was -0.39. Furthermore, the correlation coefficients for T cells and M1 macrophages were 0.48, 0.55, and 0.35, respectively, and the correlation coefficient for M0 macrophages was -0.43. The ultra-high positive correlation between immature B cells and plasma cells reached 0.78 (Figure 5C).

Figure 5. Cibersort analysis results. (A) The relative proportion of immune cells in 84 osteosarcoma patients; (B) The violin diagram shows the difference in the proportion of immune cell infiltration between the low-risk group and the high-risk group. (C) Correlation between immune cells.

According to the relative proportion of 22 immune cells, TCGA osteosarcoma patients were classified into high- and low-risk groups. The impact of various types of immune cells on the survival of osteosarcoma was examined by employing Kaplan-Meier survival analysis. The outcomes highlighted that individuals with a higher proportion of CD8+T cells and activated CD4 memory T cells experienced a more favorable prognosis. On the other hand, patients with a higher proportion of resting dendritic cells, activated mast cells, and immature CD4 T cells demonstrated a poorer prognosis (Figure 6). In addition, the link between oxidative phosphorylation-related differentially expressed genes (DEGs) and immune cells in risk scores and five risk-scoring models was compared. The risk score exhibited a positive link to activated CD4 naive T cells and mast cells while showing an inverse link to activated CD4 memory T cells, CD8+T cells, and follicular helper T cells (Supplementary Figure 5). The expression of the ATP6V0D1 gene was proportional to M1 macrophages, CD8+T cells, and M2 macrophages, and inversely proportional to activated mast cells, M0 macrophages, and CD4 naive T cells (Supplementary Figure 6). Moreover, LHPP gene expression was positively correlated with CD8+T cells and follicular helper T cells, while negatively correlated with CD4 naive T cells (Supplementary Figure 7). The expression of the COX6A2 gene was proportional to naive B cells and inversely proportional to memory B cells (Supplementary Figure 8). MTHFD2 gene expression was not linked to the 22 immune cells examined in this study. NDUFB9 gene expression was proportional to activated mast cells and monocytes (Supplementary Figure 9). Furthermore, the low-risk group exhibited statistically higher matrix, immune, and ESTIMATE scores (matrix score plus immune score). Moreover, the low-risk group exhibited a higher PDCD1LG2 expression, suggesting that this group might be more responsive to PD1 combined with immunosuppressive therapy (Supplementary Figure 10).

Figure 6. Survival analysis of immune cells. The impact of Dendritic cells resting (A), Mast cells activated (B), T cells CD4 memory activated (C), T cells CD4 naïve (D) and T cells CD8 (E) on the survival of osteosarcoma were examined by employing Kaplan-Meier survival analysis.

Drugs with potential therapeutic effects on osteosarcoma

In order to further explore why oxidative phosphorylation-related genes in osteosarcoma were significantly correlated with the prognosis of osteosarcoma, this study conducted a drug sensitivity analysis of osteosarcoma. Drug sensitivity analysis showed that Bortezomib (p = 1.3e − 06), CGP. 082996 (p = 7.7e − 05), GNF. 2 (p = 0.00013), MG.132 (p = 3.2e − 06), NVP.TAE684 (p = 1.8e − 05), PAC.1 (p = 0.00018), PF.02341066 (coxotinib) (p = 3.4e − 06), and Roscovitine (p = 2.5e − 05) demonstrated substantial sensitivity in both risk groups. Individuals in the low-risk group displayed more sensitivity to Bortezomib, CGP. 082996, GNF. 2, MG. 132, NVP. TAE684, F.02341066, and Roscovitine. In contrast, the individuals at higher risk exhibited greater sensitivity to PAC (Supplementary Figure 11). The stability of proteins that hinder cell survival and cell cycle advancement, such as p53, was increased by bortezomib. CGP60474 is a protein kinase C (PKC) inhibitor, and GNF-2 is a Bcr-Abl fusion gene inhibitor. MG-132 is an inverse proteasome inhibitor with the ability to induce tumor cell apoptosis. Moreover, NVP-TAE684 is a selective inhibitor of anaplastic lymphoma kinase (ALK), which is associated with the pathogenesis of various cancers and can serve as an important therapeutic target. PAC-1 drugs have the ability to induce cell apoptosis, and their use alone or in combination with chemotherapy has shown anticancer effects in lung cancer, melanoma, and osteosarcoma, among others. Furthermore, PF-02341066 was an ALK inhibitor that can induce autophagy in various tumor cell lines by inhibiting the STAT3 pathway. Roscovitine was a selective inhibitor of cyclin-dependent kinase (CDK). CDK is an essential cell cycle regulator and often participates in the deregulation of human tumors.

Nomogram of the prediction model

Univariate and multivariate Cox regression analyses were conducted to explore the link between gender, age, metastasis, risk-scoring, and prognosis of individuals who suffer from osteosarcoma (Supplementary Figure 12A, 12B). The findings demonstrated that both metastasis and risk score were identified as independent predictors of prognoses for individuals suffering from osteosarcoma. The 1-, 3-, and 5-year OS of the patients were predicted using a nomogram that took the outcomes of multivariate Cox regression into account (Figure 7A). This nomogram included the risk-scoring index of the patient and whether the patient had metastasis. The predicted performance results of the nomogram estimated by the ROC curve highlighted that the respective AUC values of 1-, 3-, and 5-year OS of the training set were 0.929, 0.807, and 0.809 (Figure 7B). However, the respective AUC values of 1-, 3- and 5-year OS of the verification set were 0.786, 0.885, and 0.924 (Figure 7C). The C index also highlighted a good predictive performance exhibited by the prognosis model (Figure 7D, 7E).

Figure 7. Validation of the prognostic model nomogram. (A) Nomogram. (B) Time ROC curve of training set. (C) Time ROC curve of the validation set. (D) The C-index calibration curve for age, gender, risk score, and nomogram of the training set. (E) Validate the C-index calibration curve for age, gender, risk score, and nomogram.

Gene features set enrichment analysis

To investigate the possible mechanism underlying the association between oxidative phosphorylation-related genes and the prognosis of individuals with osteosarcoma, an enrichment analysis of DEGs (|log2FC|>1, p-value<0.05) was carried out across the two risk groups. GO enrichment analysis revealed that the DEGs were primarily associated with biological processes related to myofibril assembly and muscle contraction. In terms of cellular components, the enrichment was predominantly observed in myofibrils. Furthermore, the molecular functions of these genes were found to be primarily associated with muscle structural components. The pathways for KEGG enrichment analysis included neural active ligand- receptor interaction pathways. GSEA enrichment analysis highlighted that the enrichment pathways of DEGs between both risk groups included leukocyte adhesion, cell activation, actin-mediated cell contraction, and positive regulation of external stimulus responses (Supplementary Figure 13).

Data source and preprocessing

Osteosarcoma transcriptome and clinical data were retrieved from the UCSC Xena (http://xena.ucsc.edu/) as the training set. A validation set for this study was obtained by acquiring transcriptome and clinical data of individuals diagnosed with osteosarcoma from the GSE21257 datasets. The datasets were retrieved from the publicly accessible GEO (Gene Expression Omnibus) database. Researchers can access these datasets by following the provided link: https://www.ncbi.nlm.nih.gov/geo. The GSE28425 dataset was used to screen for osteosarcoma differential genes.

Development of the risk-scoring model

Using the oxidative phosphorylation-related genes screened out above, the model was constructed based on the training set data. Differential genes above 0.5-fold were screened using GSE28425, whereas univariate Cox regression analysis was conducted to screen prognostic oxidative phosphorylation-related genes. Subsequently, LASSO regression analysis was conducted to calculate the risk coefficient (coefi) of all oxidative phosphorylation-related genes. To compute the risk score for every individual, the expression (expi) and coefficient (coefi) values of all oxidative phosphorylation genes related to prognosis were added. This was done for every sample in the study, and the resulting value was considered as the risk score of the patient (risk score = (coefi × expri)). By utilizing the median value of risk scores, individuals were classified into high- and low-risk groups. In order to examine whether a noteworthy disparity in the overall survival of individuals between both groups existed, a Kaplan-Meier survival analysis was conducted. The capability of the risk scoring model for making predictions was analyzed by calculating the area under the curve (AUC) using receiver operating characteristic (ROC) curve analysis.

Analysis of subgroups and confirmation of the capacity of the risk prognostic model for independent prognostication

Patient subgroups were created based on age (<18, ≥ 18), sex (male, female), and if metastasis had taken place (metastatic or non-metastatic). To ascertain whether there were appreciable variations in OS across high- and low-risk patients in each subgroup, Kaplan-Meier survival analysis was utilized.

Analysis of tumor microenvironment across the two risk groups

The abundance of 22 types of immune cells in both risk groups (high and low) was calculated by utilizing the CIBERSORT algorithm (https://cibersortx.stanford.edu/). Additionally, high and low groups were established as per the median expression of different immune cells, and the difference in survival time across the two groups was determined. The ESTIMATE (Estimation of Stromal and Immunological Cells in Malignant Tumor Tissues using Expression Data) algorithm was employed to assess the stromal, immunological, and ESTIMA scores, as well as the tumor purity of various molecular subgroups. Ultimately, the association of genes with immune cells was examined.

Drug sensitivity analysis

Drug sensitivity analysis was conducted using “Rs” “limma”, “ggpubr”, and “pRRophic” packages to examine the varying sensitivities of drugs for patients in both risk groups in the bone and meat data of TCGA. Moreover, the potential osteosarcoma treatment drugs were screened with a screening criterion of p-value < 0.001.

Nomogram of the prediction model

According to the clinical information of sex, age, metastasis, and risk score of the oxidative phosphorylation gene, univariate and multivariate Cox regressions were implemented, and the nomogram was constructed according to the results. Following the calculation of OS for all patients, their 1-year, 3-year, and 5-year survival rates were predicted, and a ROC curve was developed to measure the predictive accuracy of the nomogram.

Gene set enrichment analysis

In addition, enrichment analyses (GSEA, GO, and KEGG) were conducted on the differentially expressed genes (DEGs) identified between both risk groups.

Statistical analysis

R (v 4.2.0) was employed to carry out all statistical analyses. OS curves were generated by utilizing the Kaplan-Meier method and assessed by means of the log-rank test. A significance level of 0.05 or greater was used for all two-tailed p-values.

Data availability statement

The data that support the findings of this study are openly available in Therapeutically Applicable Research to Generate Effective (TARGET) database and GEO (Gene Expression Omnibus).