Probing the Potential of Defense Response-Associated Genes for Predicting the Progression, Prognosis, and Immune Microenvironment of Osteosarcoma

Abstract

Background: The defense response is a type of self-protective response of the body that protects it from damage by pathogenic factors. Although these reactions make important contributions to the occurrence and development of tumors, the role they play in osteosarcoma (OS), particularly in the immune microenvironment, remains unpredictable.

Methods: This study included the clinical information and transcriptomic data of 84 osteosarcoma samples and the microarray data of 12 mesenchymal stem cell samples and 84 osteosarcoma samples. We obtained 129 differentially expressed genes related to the defense response (DRGs) by taking the intersection of differentially expressed genes with genes involved in the defense response pathway, and prognostic genes were screened using univariate Cox regression. Least absolute shrinkage and selection operator (LASSO) penalized Cox regression and multivariate Cox regression were then used to establish a DRG prognostic signature (DGPS) via the stepwise method. DGPS performance was examined using independent prognostic analysis, survival curves, and receiver operating characteristic (ROC) curves. In addition, the molecular and immune mechanisms of adverse prognosis in high-risk populations identified by DGPS were elucidated. The results were well verified by experiments.

Result: BNIP3, PTGIS, and ZYX were identified as the most important DRGs for OS progression (hazard ratios of 2.044, 1.485, and 0.189, respectively). DGPS demonstrated outstanding performance in the prediction of OS prognosis (area under the curve (AUC) values of 0.842 and 0.787 in the training and test sets, respectively, adj-p < 0.05 in the survival curve). DGPS also performed better than a recent clinical prognostic approach with an AUC value of only 0.674 [metastasis], which was certified in the subsequent experimental results. These three genes regulate several key biological processes, including immune receptor activity and T cell activation, and they also reduce the infiltration of some immune cells, such as B cells, CD8+ T cells, and macrophages. Encouragingly, we found that DGPS was associated with sensitivity to chemotherapeutic drugs including JNK Inhibitor VIII, TGX221, MP470, and SB52334. Finally, we verified the effect of BNIP3 on apoptosis, proliferation, and migration of osteosarcoma cells through experiments.

Conclusions: This study elucidated the role and mechanism of BNIP3, PTGIS, and ZYX in OS progression and was well verified by the experimental results, enabling reliable prognostic means and treatment strategies to be proposed for OS patients.

Keywords: defense response; immune; metastasis; osteosarcoma; prognosis; therapy.

Figure 1

Screening of DEGs related to the defense response. (A) Volcano plot illustrating the DEGs in osteosarcoma and normal groups with the threshold set at |logFC| ≥ 1 and adj-p ≤ 0.05. (B) DEGs are significantly enriched in the GOBP_DEFENSE_RESPONSE pathway. (C) Trend in the number of studies on GOBP_DEFENSE_RESPONSE pathways in recent years. (D) DRGs obtained by taking the intersection of DEGs and GOBP_DEFENSE_RESPONSE pathway genes. (E) Heatmap showing the expression of DRGs in osteosarcoma samples and normal samples.

Figure 2

Obtaining DRGs associated with osteosarcoma prognosis. (A) Univariate Cox regression analysis for identifying prognostic DRGs. (B) Interaction network diagram of prognosis-related DRGs. (C,D) Lasso–Cox regression analysis was performed to construct prognostic prediction models.

Figure 3

Evaluation of DGPS. (A) Univariate Cox analysis. Risk score and metastasis were statistically significant. (B) Multivariate Cox analysis. (C) ROC curve of DGPS in training group. (D) ROC demonstrating that the predictive accuracy of DGPS was superior to the other clinical parameters in the training set. (E) Kaplan–Meier curves of overall survival in the training set. (F,G) Distribution of risk scores and distribution of overall survival status and risk score in the training set. Blue: low risk; red: high risk. (H) Heatmap indicating the expression degrees of BNIP3, PTGIS, and ZYX in the training set. ROC curve, receiver operating characteristics curve; AUC, area under the curve; p < 0.05, statistically significant.

Figure 4

Kaplan–Meier plots depicting subgroup survival analyses stratified by gender, age, and metastasis.

Figure 5

Verification of DGPS. ROC curve of DGPS in the test set (A) and in the entire cohort (B). ROC demonstrated that the predictive accuracy of DGPS was superior to that of other clinical characteristics in the test set (C) and in the entire cohort (D). Kaplan–Meier curves of overall survival (OS) in the test set (E) and in the entire cohort (F). Survival status of patients with osteosarcoma in the test set (G,I) and in the entire cohort (H,J). Blue: low risk; red: high risk. The heatmap indicates the expression degrees of BNIP3, PTGIS, and ZYX in the test set (K) and in the entire cohort (L).

Figure 6

PCA plots depicting the distribution of samples based on the expression of model genes (A), DRGs (B), and all genes (C). Differential expression of model genes in the high- and low-risk groups is shown in box plots (D–F).

Figure 7

Construction and evaluation of a nomogram based on DGPS. Nomogram used to predict prognosis was constructed based on DGPS in the training set (A), test set (E), and entire cohort (I). Calibration curves of the nomogram in the training set (B), test set (F), and entire cohort (J). The C-index curves for assessing the discrimination ability of DGPS and other clinical characteristics at each time point in the training set (C), test set (G), and entire cohort (K). ROC curves of the nomograms at one, three, and five years in the training set (D), test set (H), and entire cohort (L). “*”represented “p < 0.05”, “***”represented “p < 0.001”.

Figure 8

Exploration of the association of tumor metastasis with BNIP3. (A) Correlations between BNIP3 and osteosarcoma metastasis are displayed in box plots. ROC curve of diagnosis of osteosarcoma metastasis by BNIP3 in the training set (B) and in the entire cohort (C). (D) Relationship between the expression of BNIP3 and the expression of tumor metastasis-related genes MYC, NELL1, SAR1A, PLOD2, TNFAIP8L1, and TRIM22. (E,F) Expression of BNIP3 in different tumors. “*”represented “p < 0.05”.

Figure 9

Pathway enrichment analysis of the genes most strongly associated with BNIP3. (A) Heatmap showing the 20 genes with the strongest positive or negative correlations with BNIP3 expression. (B,C) Pathway enrichment analysis showing the enrichment of genes in different pathways.

Figure 10

GO and KEGG pathway enrichment analyses. (A) Bar plot of the top 10 GO enrichment terms. (B) Bar plot of the top 30 KEGG enrichment terms. (C) Bubble chart of the top 10 GO enrichment terms. (D) Bubble chart of the top 30 KEGG enrichment terms. (E) Circle diagram of GO enrichment analysis. (F) Circle diagram of KEGG enrichment analysis. GO enrichment terms include biological process, cellular component, and molecular function.

Figure 11

Immunoassay showing that DGPS is closely related to the immune system. (A) Analysis of TMB differences between high- and low-risk groups of patients with osteosarcoma. Box plots of the ssGSEA scores of 15 immune checkpoints (B), 13 immune cells (C), and 12 immune-related functions (D) between different risk groups. (E) Heatmap showing the landscape of immune characteristics and the tumor microenvironment in the TARGET cohort determined by the ssGSEA algorithm. “*”represented “p < 0.05”,“**”represented “p < 0.01”,“***”represented “p < 0.001”.

Figure 12

The association between immune functions and risk scores and immune function scores between different risk subgroups.

Figure 13

The association between immune cells and risk scores and immune cell scores between different risk subgroups.

Figure 14

Correlation analysis of immune-related scores and risk scores. (A–C) Analysis of the variability of risk scores among different StromalScore (A), ImmuneScore (B), and ESTIMATEScore (C) subgroups. (D–F) Scatter plots of correlations between risk scores and StromalScore (D), ImmuneScore (E), and ESTIMATEScore (F).

Figure 15

Drug correlation and sensitivity analyses with JNK Inhibitor VIII, TGX221, MP470, and SB52334.

Figure 16

BNIP3 regulates the apoptosis of osteosarcoma cells. (A–E) Apoptosis of osteosarcoma cells after knockdown or overexpression of BNIP3. (F) Overexpression of BNIP3 inhibits apoptosis of osteosarcoma cells, while knockdown of BNIP3 promotes apoptosis of osteosarcoma cells. “NS” represented “No significant difference”, “***” represented “p < 0.001”.

Figure 17

BNIP3 regulates the proliferation of osteosarcoma cells. Knockdown of BNIP3 inhibits the proliferation of osteosarcoma cells, while overexpression of BNIP3 promotes the proliferation of osteosarcoma cells. “NS” represented “No significant difference”, “*” represented “p < 0.05”, “**” represented “p < 0.01”, “***” represented “p < 0.001”.

Figure 18

BNIP3 regulates the migration ability of osteosarcoma cells; knockdown of BNIP3 inhibits their migration ability, while overexpression of BNIP3 promotes the migration ability of osteosarcoma cells. “NS” represented “No significant difference”, “**” represented “p < 0.01”, “***” represented “p < 0.001”.