Identification of Gene as Predictive Biomarkers for the Occurrence and Recurrence of Osteosarcoma

Abstract

Purpose

Osteosarcoma is the most common malignant bone cancer affecting adolescents and young adults. This study aimed to screen potential diagnostic and therapeutic markers for osteosarcoma.

Methods

Differential expression analysis between osteosarcoma and control was performed in GSE99671, the differentially expressed genes (DEGs) were subjected to co-expression analysis. Enrichment analysis was employed to identify the biological functions and KEGG signaling pathways of module genes. In addition, a differential analysis was also performed between recurrent and non-recurrent osteosarcoma samples in GSE39055, and enrichment analysis was performed for DEGs. Further, Kaplan–Meier curve analysis was performed on the module genes, and receiver operating characteristic (ROC) curve was drawn. Comparison of the module with the highest correlation to osteosarcoma identified key genes. Cox regression model was utilized to identify the predictive ability of key genes for the prognosis of osteosarcoma.

Results

A total of 13 co-expression modules were identified from 4871 DEGs of GSE99671, module 1 had the highest positive correlation with osteosarcoma. Module genes were mainly enriched in autophagy and macrophage migration functions. A total of 1126 DEGs were obtained from GSE39055, significantly involved in neutrophil mediated immunity. Screening of genes with area under the ROC curve (AUC) values greater than 0.73 in both GSE99671 and GSE39055 identified 5 key genes when compared with genes from module 1. The nomogram results showed that ATF5, CHCHD8, ENOPH1, and LOC286367 might predict 5-year or 8-year survival time of osteosarcoma patients. The Cox model results confirmed that the signals of ATF5, CHCHD8, and LOC286367 were robust, and it may be used in the diagnosis, treatment, and prognosis of osteosarcoma.

Conclusion

We found that ATF5, CHCHD8, and LOC286367 can effectively identify osteosarcoma tumorigenesis and even recurrence status. This is helpful for early diagnosis and treatment, improving the clinical treatment of patients with osteosarcoma.

Keywords:

• osteosarcoma
• predictive biomarkers
• recurrence
• weighted co-expressed network analysis

Introduction

Osteosarcoma is the most common primary bone malignancy, mainly affecting children and young adults.¹ Osteosarcoma consists of malignant osteoblasts producing immature bone and bone tissue.² Although standard treatment with surgical resection and adjuvant chemotherapy has significantly improved the 5-year survival rate of osteosarcoma patients to approximately 60–70%, no significant progress has been made in improving the survival rate of patients with recurrence or metastasis over the past 30 years.^3,4 The lack of understanding of the molecular mechanisms underlying the occurrence and recurrence of osteosarcoma has severely hampered improved patient survival. When diagnosed, 40% of metastases occur in patients with advanced osteosarcoma.⁵ Therefore, elucidating the functions of osteosarcoma-related genes and exploring the possible pathological mechanisms of osteosarcoma initiation, development and recurrence are crucial for the future detection and treatment of osteosarcoma.

Depending on the histological morphology, three main categories can be distinguished: high-grade, which includes most subtypes, and intermediate and low-grade, which include periosteal and periosteal.⁶ Conventional osteosarcoma refers to high-grade tumors with intramedullary growth and is the most common type, accounting for 85% of all osteosarcoma cases during childhood and adolescence.⁷ The osteosarcoma tumor microenvironment is composed of osteosarcoma cells, osteocytes, stromal cells, vascular cells, immune cells, and the extracellular matrix (ECM).⁸ This creates a complex environment for tumor growth. More immune infiltration is found in osteosarcoma tumors to promote a local immune tolerant environment.⁹ Given the close connection between bone tissue and the immune system, it has been speculated that osteosarcoma may use similar mechanisms to evade immune recognition.¹⁰

With the development of molecular biology technology, tumor gene therapy for osteosarcoma has potential clinical applications.^11,12 Accumulating evidence indicates that the occurrence, development and prognosis of osteosarcoma are closely related to molecular mechanisms.^13,14 Through high-throughput sequencing, gene expression in osteosarcoma can be compared to normal samples, leading to an initial selection of potential targets for anticancer therapy.¹⁵ This also includes some common long non-coding RNAs (lncRNAs).¹⁶

This study aimed to discover the occurrence- or recurrence-related potential markers according to the osteosarcoma-related gene expression profiles in Gene Expression Omnibus (GEO) database. Multiple differentially expressed genes (DEGs) were screened using the weighted gene co-expression network analysis (WGCNA) algorithm. The main signaling pathways of osteosarcoma were analyzed by Gene Ontology (GO), function and Kyoto Encyclopedia of genes (KEGG) pathways. Subsequently, Kaplan–Meier and Cox models were utilized to screen the key genes related to the prognosis of osteosarcoma. The findings may provide new biomarkers and therapeutic target molecules for the occurrence and recurrence of osteosarcoma.

Materials and Methods

Data Sources

Osteosarcoma data were collected from the gene expression omnibus (GEO) (https://www.ncbi.nlm.nih.gov/geo) databases. GSE99671 included mRNA expression profiling of fresh-frozen bone samples from 18 tumoral samples and 18 non-tumoral paired samples by high throughput sequencing.^17,18 The raw data were standardized and normalized through R package DEseq2;¹⁹ then, the gene expression profile was provided in Table S1. GSE39055 included mRNA expression profiling of formalin-fixed, paraffin-embedded (FFPE) samples from 37 unique diagnostic biopsy specimens of osteosarcoma with (n=18) or non-recurrence (n=19) by array.²⁰ Expression data were processed using the R package lumi²¹ and provided in Table S2. All data were obtained from an open-access database, thus, acquiring ethical approval was not necessary.

Differential Analysis

Next, the DEseq2 package was used to identify differentially expressed genes (DEGs) in osteosarcoma (n = 18) and healthy controls (n = 18) in GSE99671 data. The limma package²² was used to screen DEGs between the recurrent osteosarcoma and non-recurrent osteosarcoma in GSE39055 data. P < 0.05 was set as the screening condition. The expression of DEGs was visualized with volcano plots.

Construction of Networks

The coexpression network for DEGs in GSE99671 was performed using Weighted correlation network analysis (WGCNA) by “WGCNA” R package.²³ The samples were used to construct scale-free topology networks, and all gene adjacencies were calculated to make a topological overlap matrix (TOM). The soft-thresholding power was chosen and used as the correlation coefficient threshold. Then, a minimum number of genes in the modules were built. The expression pattern of eigengene in each module is condensed into “module eigengene (ME)”. Genes in MEs were considered had similar expression patterns. The module–trail relationships were demonstrated using Pearson correlation analysis.

Enrichment Analysis

The biological process (BP) in Gene Ontology (GO) and Kyoto Encyclopedia of Genes and Genomes (KEGG) for module genes or DEGs in GSE39055 were performed using R software package clusterProfiler.²⁴ Gene set variation analysis (GSVA) was carried out using the GSVA package.²⁵ For each sample, a score for the enrichment of a set of genes using gene expression profile was obtained. A P value < 0.05 was considered statistical significance. The R package clusterProfiler was used to obtain the background set for gene set enrichment analysis (GSEA). GSEA runs in Java environment and conducted between osteosarcoma and control samples.

Screening of Survival Related Genes

Kaplan–Meier (K-M) survival analysis with the R package survival was applied to identify overall survival (OS) associated genes among the module genes. Prognostic relevant genes among the key genes were determined using Cox regression analysis in a survival package based on gene expression values and survival status data. P < 0.05 was considered to be a statistically significant difference. Based on the correlation coefficients of the risk model, a formula for calculating the prognostic risk score for each patient was established. According to the median value of the risk score, patients were divided into high- and low-risk groups.

ROC Curve

To evaluate whether the identified survival-related genes have significant diagnostic value for osteosarcoma or recurrent osteosarcoma, we performed ROC analysis using the R software package pROC.²⁶ The area under the ROC curve (AUC) was calculated for each gene. The diagnostic accuracy of key biomarkers was evaluated using AUC values. Here, this gene can distinguish osteosarcoma from healthy individuals or recurrent osteosarcoma from non-recurrent osteosarcoma when the AUC value is greater than 0.6.

Single Sample Gene Set Enrichment Analysis (ssGSEA)

The marker gene expression information for immune cell types was obtained by Bindea et al.²⁷ The infiltration scores of immune cells were calculated using ssGSEA in GSVA R software package. Correlations between immune cells and key genes were calculated using Pearson correlation. P value <0.05 was considered significant.

Results

Coexpression of Differentially Expressed Genes

The flowchart of this study is shown in Figure 1. To identify differences in gene expression between osteosarcoma and controls, we performed a differential analysis of the osteosarcoma and control groups in GSE99671. Using threshold screening, we obtained 4871 differentially expressed genes (DEGs) between osteosarcoma and controls (Figure 2A). The co-expression behaviors of these genes were identified by performing WGCNA on the DEGs. To ensure that the coexpression network could obey the scale-free criterion, we selected β = 9 as a soft-thresholding (Figure 2B). We identified a total of 13 coexpression modules (Figure 2C). Pearson correlation analysis showed that MEturquoise (M1) had the highest positive correlation with osteosarcoma (Figure 2D).

Figure 1 The flowchart of this study. The GSE99671 data and GSE39055 data were used to identify potential biomarkers related to the occurrence and recurrence of osteosarcoma. The module genes identified in GSE99671 were used to evaluate their prognostic and diagnostic value in GSE39055. Further screening of key genes may be a target for diagnosis and treatment of osteosarcoma.

Abbreviations: WGCNA, weighted gene co-expression network analysis; GSEA, gene set enrichment analysis; ROC, receiver operating characteristic curve; AUC, area under the ROC curve.

Figure 2 Coexpression modules of differentially expressed genes. (A) The differentially expressed genes between osteosarcoma and controls in GSE99671. (B) The scale-free fit index and the mean connectivity for various soft-thresholding powers. (C) Common genes were identified for thirteen coexpression modules by WGCNA. (D) Correlations between modules and clinical trait were analyzed by Pearson correlation.

Biological Functions of Module Genes

To identify the biological roles of module genes, we performed enrichment analysis. It was found that module genes were significantly enriched in activated “positive regulation of macrophage migration“, “chaperone mediated autophagy“, “MHC protein complex assembly“ biological processes; and inhibited “neutrophil aggregation“, “regulation of vascular permeability“, “regulation of cell killing“ in osteosarcoma (Figure 3A). Using KEGG enrichment results, we found that activated HIF-1 signaling pathway, PI3K-Akt signaling pathway, and autophagy animal were significantly enriched by module genes, and inhibited B cell receptor signaling pathway, Rap1 signaling pathway, and ECM receptor interaction were also significantly enriched in osteosarcoma (Figure 3B). GSEA results showed that protein export, RNA polymerase, and proteasome were significantly enriched in osteosarcoma, primary immunodeficiency, PPAR signaling pathway, and neutrophil extracellular trap formation were significantly enriched in control samples (Figure 3C).

Figure 3 Biological functions of module gene enrichment. (A) Module genes significantly enriched for activated and inhibited biological processes. (B) Module genes significantly enriched for activated and inhibited KEGG pathway. (C) The activated and inhibited KEGG pathways in GSEA of module genes involved. P < 0.05 was considered statistically significant.

Abbreviation: GSVA, gene set variation analysis.

Recurrence Related Changes of Osteosarcoma

To identify molecular alterations associated with osteosarcoma recurrence, we performed a differential analysis of recurrent and non-recurrent osteosarcoma gene expression in GSE39055. A total of 1126 differentially expressed genes were obtained (Figure 4A). In GO enrichment results, the activated “neutrophil mediated immunity“, “skeletal muscle fiber differentiation“, “response to iron (III) ion“ and inhibited “neutrophil apoptotic process“, “leukocyte homeostasis“, “mitochondrial organization“ were significantly enriched by differentially expressed genes (Figure 4B). KEGG enrichment results identified that the differentially expressed genes were significantly involved in the activated “endocrine resistance“, “oxidative phosphorylation“, “renin angiotensin system“; and inhibited “p53 signaling pathway“, “cell cycle“, “FOXO signaling pathway“ in recurrent osteosarcoma (Figure 4C).

Figure 4 The potential molecular changes in osteosarcoma recurrence. (A) Volcano plot of differentially expressed genes between recurrent osteosarcoma and non-recurrent osteosarcoma. Red for activated and green for inhibited. The four genes with the largest fold change of activated or inhibited were labeled. (B) Differentially expressed genes involved in activated and inhibited biological processes quantified by gene set variation analysis (GSVA). The longer the column, the more genes involved in this term. (C) Differentially expressed genes involved in activated and inhibited KEGG pathway quantified by gene set variation analysis (GSVA). The longer the column, the more genes involved in this term.

Abbreviation: BP, biological processes.

Identification of Potential Diagnostic Genes

From the survival information of 37 osteosarcoma patients obtained from GSE39055, we found that the patient’s 5-year overall survival (OS) was 35% (Figure 5A). The survival probability increased per year already survival related to the total survival time. To identify key genes involved in the prognosis of osteosarcoma, we performed KM curve analysis of module genes. Among the module genes, we identified 262 genes that significantly affected osteosarcoma survival. By plotting the ROC curves of these survival-related genes, we screened for genes with AUC values greater than 0.73 in both GSE99671 and GSE39055 (Figure 5B). Including AEBP1, ATF5, CHCHD8, DYRK3, ENOPH1, GMIP, LOC286367, PKP4, R3HDM1, TRIM66, and ZMYND17. They were also differentially expressed genes in GSE39055 and thus may predict both osteosarcoma occurrence and recurrence. In addition, five genes (ATF5, CHCHD8, ENOPH1, LOC286367, and R3HDM1) were identified as potentially most strongly associated with osteosarcoma by contrasting with the gene for module 1. They were identified as key genes. Nomograms were constructed using Cox regression analysis results, and differential expression of ATF5, CHCHD8, ENOPH1, and LOC286367 could predict 5 - and 8-year OS of osteosarcoma patients (Figure 5C).

Figure 5 Screening of key genes predicting the occurrence and recurrence of osteosarcoma. (A) Kaplan–Meier estimates for conditional survival up to 6 years in 37 patients given 0–5 years’ survival of osteosarcoma. Each column represents the survival time, and each row represents the percentage reaching the specified survival time. (B) On the left are survival associated genes with AUC values in GSE99671 and GSE39055 and on the right are genes with AUC values greater than 0.73 in both datasets. Red represents up-regulated expression and blue represents down regulated expression. The length of the column represents the mean AUC value. (C) Nomogram for the prediction of overall survival to achieve 5-year or 8-year survival time.

Abbreviation: AUC, area under the ROC curve.

On the other hand, osteosarcoma patients in GSE39055 were divided into high-risk and low-risk groups according to the median risk score which calculated by the Cox model (Figure 6A). The expressions of key genes in the high-risk and low-risk groups were different. Patients with death status were significantly higher in the high-risk group than in the low-risk group. In addition, we calculated the correlation of key genes and immune cells, LOC286367 was positively correlated with most immune cells, other key genes were negatively correlated with more immune cells (Figure 6B). The prediction capability of the key genes was evaluated by calculating the area under the ROC curve (AUC). The AUC values of ATF5, CHCHD8, and LOC286367 for predicting OS were greater than 0.6 in the first, third, and fifth year of osteosarcoma, indicating that they had good performance (Figure 6C). High expression of ATF5, CHCHD8, and LOC286367 was associated with the worst OS in osteosarcoma patients (Figure 6D). ATF5, CHCHD8, and LOC286367 were expressed at significantly higher levels of osteosarcoma compared with the normal group in GSE99671 (Figure 6E).

Figure 6 Diagnostic role of key genes predicting the survival of osteosarcoma. (A) Distribution of risk score, overall survival, and overall survival status and heatmap of the five key genes in the GSE39055. (B) Correlations between immune infiltrating cells and key genes were calculated using Pearson. *P < 0.05, **P < 0.01. (C) Time-dependent ROC curves measuring the predictive value of key genes in GSE39055 for 3-year, 5-year or 10-year survival time. (D) Effect of key genes expression on overall survival by Kaplan–Meier analysis in 37 patients with osteosarcoma. (E) The expression levels of ATF5, CHCHD8, and LOC286367 in osteosarcoma and normal of GSE99671. ***P < 0.001.

Abbreviation: OS, overall survival.

Discussion

Osteosarcoma is one of the most common aggressive bone tumors and is currently treated with chemical drugs combined with surgical resection. A major unsolved problem is the poor prognosis characterized by drug resistance, recurrence and metastasis.²⁸ The identification of gene signatures is crucial both for a better understanding of the molecular basis of osteosarcoma progression and for the discovery of novel targets.²⁹ In the present study, we focused on the potential target genes with prognostic value for osteosarcoma. Gene expression profiling of differentially expressed genes detected by GSE99671 in osteosarcoma tissue samples established 13 coexpression modules. Module genes were mainly enriched in immune, inflammatory and metabolic responses. In the gene co-expression network, module 1 was most significantly associated with osteosarcoma, in which genes significantly affecting survival may be potential target genes. Furthermore, combining GSE39055, we identified 5 key genes that may serve as diagnostic markers for the occurrence and recurrence of osteosarcoma. Finally, high expression of ATF5, CHCHD8, and LOC286367 was associated with worse prognosis in osteosarcoma.

WGCNA is a systems biology approach that describes correlation patterns between genes in transcriptome samples with soft threshold algorithms.²³ The results of GO and KEGG pathway enrichment analysis of the module genes led us to focus on the biological functions of autophagy and macrophage migration, as well as the HIF-1 signaling pathway and PI3K-Akt signaling pathway. Autophagy promotes the proliferation and development of osteosarcoma cells and resists tumor treatment.³⁰ Autophagy may be involved in drug sensitivity or chemoresistance during osteosarcoma treatment.³¹ Macrophages are an important immune component in the osteosarcoma microenvironment. Macrophages are highly plastic and the inflammatory phenotype (M1) and anti-inflammatory phenotype (M2) may play opposite roles in the progression of osteosarcoma.³² Activation of the HIF-1 signaling pathway promotes osteosarcoma cell growth and is a promising therapeutic target.³³ Accumulating evidence suggests that the PI3K/Akt pathway is involved in cancer initiation and progression, such as tumorigenesis, apoptosis inhibition, proliferation and drug resistance.³⁴

To further identify the underlying molecular mechanisms of osteosarcoma recurrence, we performed enrichment analysis of the differentially expressed genes between recurrence and non-recurrence. We found multiple immune related pathways, neutrophil mediated immunity, neutrophil apoptotic process, and leukocyte homeostasis. They may be associated with metastasis and recurrence of osteosarcoma.³⁵ In addition to a large number of aberrant biological functions, FOXO could control the expression of genes involved in cell death and cell cycle arrest, exerting tumor suppressor activity.³⁶ Tumor suppressor p53 tumor cells have been reported to exert anticancer effects by inducing cell cycle arrest and apoptosis.³⁷

Of the 13 coexpression modules we identified, module 1 was found to be strongly associated with osteosarcoma. Among the module genes identified to be significantly associated with OS of osteosarcoma, 11 genes were screened with high specificity and sensitivity as potential molecular markers for predicting the occurrence and recurrence of osteosarcoma. Among them, 5 genes were genes in module 1 and were considered as key genes. In the Cox regression model, high expression of CHCHD8 or ENOPH1 may benefit the prognosis of osteosarcoma patients, and low expression of ATF5 or LOC286367 may prolong patient survival. However, the K-M curve results differed from the Cox analysis results in that patients with low expression of the CHCHD8 gene had better overall survival. The risk assessment model constructed by 5 key genes clearly distinguished the status of osteosarcoma survival and death. Based on the time-dependent ROC curves results, we identified ATF5, CHCHD8, and LOC286367 as potential diagnostic and prognostic biomarkers for clinical outcome prediction. ATF5 is considered an anti-apoptotic factor because it regulates the expression of the anti-apoptotic components BCL2 and MCL1.³⁸ Studies have reported ATF5 to promote tumor growth and survival and have been reported to be associated with recurrence of osteosarcoma.^39,40 CHCHD8 has been reported to be associated with drug resistance in gastric cancer, but the association with osteosarcoma has not been reported.⁴¹ LOC286367 was reported to be associated with inflammatory response,^42,43 the results of our analysis argued that it may be a potential marker for the occurrence and recurrence of osteosarcoma.

However, the present study has certain limitations. The first and most important is the lack of experimental validation. Second, the lack of detailed clinical data, such as chemotherapy regimens and tumor stages, limits in-depth association analyses. Finally, whether the biomarkers we identified can be applied in the clinic also needs a large number of samples to validate, which will be the focus of our future studies.

Conclusion

The present study identified potential diagnostic markers and useful therapeutic targets for osteosarcoma patients, and they may be able to predict patient prognosis. The mechanism of recurrence of osteosarcoma may be associated with neutrophil immunity and cell cycle arrest. In-depth exploration of the potential target genes and molecular deregulation mechanisms to develop corresponding prevention and treatment countermeasures can achieve breakthrough progress in the prognosis of osteosarcoma.

Data Sharing Statement

Data were downloaded from the Gene Expression Omnibus (GEO) database (https://www.ncbi.nlm.nih.gov/geo/).

Acknowledgments

This study was supported by the Scientific Research Project of Guangxi Health Commission (Z20180799).

Disclosure

The authors report no conflicts of interest related to this work.