Construction and validation of a novel apoptosis-associated prognostic signature related to osteosarcoma metastasis and immune infiltration.

BACKGROUND: Apoptosis played vital roles in the formation and progression of osteosarcoma. However, no studies elucidated the prognostic relationships between apoptosis-associated genes (AAGs) and osteosarcoma.
METHODS: The differentially expressed genes associated with osteosarcoma metastasis and apoptosis were identified from GEO and MSigDB databases. The apoptosis-associated prognostic signature was established through univariate and multivariate cox regression analyses. The Kaplan-Meier (KM) survival curve, ROC curve and nomogram were constructed to investigate the predictive value of this signature. CIBERSORT algorithm and ssGSEA were used to explore the relationships between immune infiltration and AAG signature. The above results were validated in another GEO dataset and the expression of AAGs was also validated in osteosarcoma patient samples by immunohistochemistry.
RESULTS: HSPB1 and IER3 were involved in AAG signature. In training and validation datasets, apoptosis-associated risk scores were negatively related to patient survival rates and the AAG signature was regarded as the independent prognostic factor. ROC and calibration curves demonstrated the signature and nomogram were reliable. GSEA revealed the signature related to immune-associated pathways. ssGSEA indicated that one immune cell and three immune functions were significantly dysregulated. The immunohistochemistry analyses of patients' samples revealed that AAGs were significantly differently expressed between metastasis and non-metastasis osteosarcomas.
CONCLUSIONS: The present study identified and validated a novel apoptosis-associated prognostic signature related to osteosarcoma metastasis. It could serve as the potential biomarker and therapeutic targets for osteosarcoma in the future.

Introduction

Osteosarcoma, one of the most aggressive diseases, accounted for one-fifth of all primary bone malignant tumors and about 2.4% of pediatric cancers [1]. It mainly occurred in the metaphysis of long bones and was accompanied by numerous characteristics such as no specific early symptoms, high tumor heterogeneity and numerous genomic instabilities [2,3]. With the progress of modern biomedicine, especially the combination of neoadjuvant chemotherapy and radical tumor resection, the 5-year survival rate of patients had significantly increased to over 70% [4]. But for patients who developed metastasis at diagnosis or after treatment, the 5-year survival rate drastically decreased to less than 20% [5]. Currently, numerous studies tried to explore and develop novel therapeutic strategies for metastatic osteosarcoma patients. However, little progress had been achieved since the molecular mechanisms of osteosarcoma metastasis were still unknown.

Programmed cell death, also widely known as apoptosis, took part in several physiological and pathological processes including the development of embryo, elimination of abnormal cells and maintenance of the whole body homeostasis [6,7]. Dysregulation of apoptosis, such as appearance or enhancement in the wrong place or time, always triggered numerous diseases including malignant tumors [8]. What's more, alterations in apoptosis pathways related to the resistance of chemotherapy, radiotherapy and targeted therapy et al., which led to the recurrence and metastasis of cancers [9]. Thus, large numbers of studies attempted to reveal the relationships between apoptosis and cancer development. This might improve the prognosis of patients.

In the last decades, the diagnosis of osteosarcoma metastasis mainly depended on radiologic methods such as X-ray, computerized tomography (CT) and positron emission tomography-computed tomography (PET-CT). However, the low sensitivity and accuracy of these detective methods lead osteosarcoma patients to lose the best diagnostic and therapeutic opportunities. Besides, the latest precision medicine also demanded specific and sensitive biomarkers to evaluate therapeutic benefits and prognosis. Nowadays, with the progression of microarray and sequencing technology, various oncogenes had been identified and some of them could act as the prognostic biomarker in cancers. A larger amount of researchers also identified that several genes could combine as the signature to predict prognosis [10,11]. In osteosarcoma, some studies had constructed prognostic signatures related to hypoxia, energy metabolism and tumor microenvironment [12], [13], [14]. Unfortunately, no apoptosis-associated signature had been established.

In this study, we screened out the apoptosis-associated differentially expressed messenger RNAs (DEmRNAs) between metastatic and non-metastatic osteosarcoma patients in Gene Expression Omnibus (GEO) dataset (GSE21257). The prognosis-related DEmRNAs were identified through univariate cox regression analysis and the signature was established through the coefficients achieved in multivariate cox regression analysis. The relationship between gene signature and patients’ prognosis was conducted through Kaplan–Meier (KM) analysis. The nomogram was also constructed and the immune infiltration of patients with different risk scores was explored. In addition, the apoptosis-associated signature was validated in another GEO dataset (GSE39055) and the DEmRNAs in signature were also validated in osteosarcoma patients’ samples by immunohistochemistry (IHC). The experimental procedure of this research was displayed in Fig. 1. We hoped this novel apoptosis-related signature could help to predict the prognosis of osteosarcoma patients and serve as the potential targets for anti-cancer treatment.

Fig. 1

The flowchart of the study. (A). The apoptosis-associated metastatic DEmRNAs were identified in GSE21257 and MSigDB database. (B). Univariate and multivariate cox regression analyses were conducted to construct the apoptosis-associated gene signature. (C). The appilication of apoptosis-associated gene signature in clinical and validation in GSE39055. (D). The immune infiltration features related to apoptosis risk. DEmRNA: differentially expressed messenger RNA. GEO: Gene Expression Omnibus.

Materials and methods

Acquisition of datasets

The gene expressive data of osteosarcoma were searched in the GEO database (http://www.ncbi.nlm.nih.gov/geo). The selection criteria of datasets were shown as follows: (1) the osteosarcoma was confirmed by pathology; (2) the follow-up data included metastatic and prognostic information; (3) the sample size of dataset was greater than 30. Finally, GSE21257 (platform GPL10295, Illumina human-6 v2.0 expression beadchip Illumina, Inc., San Diego, CA, United States) was selected as the training group and GSE39055 (platform GPL14951, Illumina HumanHT-12 WG-DASL V4.0 R2 expression beadchip Illumina, Inc., San Diego, CA, United States) was chosen as the validation group. The clinical characteristics of patients in two datasets were shown in Supplementary Table 1. The apoptosis-related genes were acquired in the hallmark gene sets in Molecular Signature Database (MSigDB, https://www.gsea-msigdb.org/gsea/msigdb/).

Screening of apoptosis-associated DEmRNAs

The mRNA expressive data and patients’ clinical information were downloaded from each dataset. Patients were divided into metastatic and non-metastatic groups according to the follow-up data. The DEmRNAs among two groups were identified by Linear Models for Microarray Analysis (limma) package in R software. P < 0.05 and |log2fold change (FC)| > 0.585 were considered as statistically significant. The heatmap of DEmRNAs was drawn by the “pheatmap” package in R software. At last, the apoptosis-associated DEmRNAs were screened out through the intersection between DEmRNAs and apoptosis-related genes.

Construction and validation of the apoptosis-associated prognostic signature

The apoptosis-associated genes related to patients’ survival were obtained from the univariate cox regression analysis. The multivariate cox stepwise regression analysis was used to establish the prognostic signature. The akaike information criterion (AIC) was employed to minimize the number of DEmRNAs which represented the signature [15]. The risk score was calculated as follows:

According to the median value of risk scores, patients were separated into high- and low-risk groups. The ‘survival’ package in R software was employed to construct the KM survival curve of apoptosis-associated signature. The predictive reliability of signature was conducted through time-related receiver operating characteristic (ROC) curves. The apoptosis-associated signature was validated in GSE39055. P< 0.05 was considered statistically significant.

Establishment of the nomogram

To stratification different risk patients, the nomogram which integrated clinical characteristics such as age, gender and metastasis was constructed [16]. The score of each clinical factor in the nomogram was estimated by multivariate cox regression analysis. Harrell's concordance index (C-index) was acquired to evaluate the calibration of this model. The nomogram was established and validated in GSE39055 again. The package ‘rms’ was used to operate the above procedures in R software.

Exploring clinical correlations and biological functions of apoptosis-associated signature

The independent prognostic role of this apoptosis-associated gene signature was verified through univariate and multivariate cox regressions. The subgroup analyses of signature and individual genes were conducted through the follow-up information. Next, the gene set enrichment analysis (GSEA) was performed to evaluate the biological functions of gene signature. The top 15 significantly enriched pathways were visualized through ‘ggplot2’ package. The validation of this signature was conducted in GSE39055. All above procedures were operated in R software. P< 0.05 was considered as statistically significant.

Evaluation of tumor immune infiltration status related to apoptosis-associated signature

The tumor immune infiltration features of osteosarcoma were estimated through CIBERSORT deconvolution algorithm just as mentioned previously [17,18]. In short, the mRNA expression in osteosarcoma was matched with the gene signature matrix from CIBERSORT platform (https://cibersortx.stanford.edu) to generate the infiltrative features of 22 immune cells in osteosarcoma. The correlations between different immune cells were performed by ‘corrplot’ package in R software.

Next, the different immune infiltration features between low- and high-risk patients were investigated through single-sample gene set enrichment analysis (ssGSEA) in GSE21257 and GSE39055 datasets. The overlapped events with the same changes were regarded as characteristic variants in immune system. P< 0.05 was considered statistically significant.

Immunohistochemistry analysis

Eight patients from Ruijin Hospital were selected for IHC assay. All patients were diagnosed as osteosarcoma by pathology. Among these patients, five developed metastases in 3 years and the others lived without metastasis.

Each osteosarcoma sample was fixed with 10% formalin, embedded in paraffin and sliced as 4 μm thickness. The sections were deparaffinized and rehydrated by xylene and ethanol. After heat-induced antigen retrieval, endogenous peroxidase activity blockage and serum sealing, the samples were incubated with rabbit anti-IER3 (ABclonal, A8566, 1:100) and anti-HSPB1 (Proteintech, 18,284–1-AP, 1:200) antibodies at 4℃ overnight. Then slices were covered with horseradish peroxidase-coupled goat anti-rabbit secondary antibody (CST, 7074, 1:200) at room temperature for 50 min and stained with diaminobenzidine. The nuclei were counterstained by hematoxylin. At last, the samples were dehydrated and mounted. All slices were visualized by the microscope (Nikon DS-U3).

Results

Acquisition of apoptosis-associated DEmRNAs related to osteosarcoma metastasis

The RNA-seq data were downloaded from GSE21257 which included 34 metastatic and 19 non-metastatic patients. A total of 217 down-regulated and 247 up-regulated DEmRNAs related to osteosarcoma metastasis were screened out (Fig. 2A). One hundred and sixty-one apoptosis-associated genes were downloaded from the hallmark gene sets in MSigDB and then intersected with DEmRNAs obtained previously (Fig. 2B). Finally, 9 apoptosis-associated DEmRNAs were identified (Table 1).

Fig. 2

Identification of apoptosis-associated DEmRNAs related to osteosarcoma metastasis. (A) DEmRNAs related to osteosarcoma metastasis in GSE21257 database. Red indicated up-regulated DEmRNAs (P〈 0.05, FC 〉 1.5) and green indicated down-regulated DEmRNAs (P > 0.05, FC < 0.67). (B) Nine apoptosis-associated DEmRNAs were identified in GSE21257 and MSigDB databases. DEmRNA, differentially expressed messenger RNA; FC, fold change; MSigDB:Molecular Signature Database.

Table 1

The apoptosis-associated DEmRNAs in GSE21257.

Gene	Log2FC	P Value
CD14	−1.29	0.000057
CD2	−0.83	0.00082
GNA15	−0.77	0.0020
HSPB1	0.59	0.0097
IER3	−0.97	0.0030
IFNGR1	−0.71	0.000093
IGFBP6	−0.70	0.021
IL1B	−0.85	0.000077
PEA15	−0.59	0.000012

DEmRNA:differentially expressed messenger RNA.

Construction of the apoptosis-associated prognostic signature

Nine DEmRNAs obtained previously were applied to construct the apoptosis-associated prognostic signature by univariate and multivariate cox regression analyses. In univariate cox regression analysis, four genes were identified as prognosis-associated genes (Fig. 3A). Subsequently, these genes were explored by multivariate cox regression analysis and two of them were screened out to establish gene signature (Table 2). The risk score formula was constructed based on the genes' regression coefficients and expression levels:

Fig. 3

Establishment of the apoptosis-associated prognostic signature in GSE21257. (A) Univariate cox regression analysis of apoptosis-associated DEmRNAs revealed four prognostic genes. (B)The distribution and median value of risk score. (C) The heatmap of HSPB1 and IER3 in high- and low-risk groups revealed that HSPB1 was positively correlated with risk score while IER3 showed the opposite trend. (D) The scatter plot of patients’ risk scores and survival status. (E) Kaplan-Meier survival plot of patients in high- and low-risk groups. (F) Time-dependent ROC curve at 1, 3, and 5 years of apoptosis-associated prognostic signature. DEmiRNA, differentially expressed microRNA; ROC, receiver operating characteristic;IER3:Immediate Early Response 3, HSPB1:Heat Shock Protein Family B Member 1.

Table 2

The multivariate cox regression analysis of GSE21257.

Gene	Coefficient	HR	HR 95Low	HR 95High	P Value
HSBP1	0.44	1.55	0.86	2.81	0.149
IER3	−0.34	0.71	0.50	1.01	0.059

HR:hazard ratio.

Next, all patients were divided into high- and low-risk groups based on the median value of risk score (Fig. 3B). The heatmap revealed that heat shock protein family B member 1(HSPB1) was increased in the high-risk group while immediate early response 3(IER3) exhibited the opposite trend (Fig. 3C). Meanwhile, the prognosis of high-risk patients was worse than low-risk group (Fig. 3D and 3E). At last, the time-dependent ROC curves were established (Fig. 3F). The area under the curve (AUC) of 1-year, 3-year and 5-year ROC plots were 0.75, 0.74 and 0.71 respectively, which indicated the robustness of this apoptosis-associated prognostic signature.

Validation of the apoptosis-associated prognostic signature

The reliability of this gene signature was verified in GSE39055. Thirty-seven osteosarcoma patients with integrated prognostic information were involved in the subsequent analysis. According to the previous risk score formula, 18 patients were separated into the low-risk group and the rest 19 were considered as high risk (Supplementary Figure 1A). The expressive tendencies of HSPB1 and IER3 were similar to GSE21257 (Supplementary Figure 1B). The risk plot (Supplementary Figure 1C) and KM survival analysis (Supplementary Figure 1D, P = 0.023) also showed that low risk scores related to better prognosis. In addition, the time-dependent ROC curve indicated that the gene signature was reliable (Supplementary Figure 1E).

Investigation of the relationship between clinical variables and apoptosis-associated prognostic signature

Univariate and multivariate cox regression analyses were applied to investigate the predictive relationships between apoptosis-associated signature and clinical variables. The risk score was identified as the independent prognostic factor in GSE21257 (HR = 2.48, 95% CI = 1.44–4.26, P = 0.0010) and GSE39055 respectively (HR = 4.69, 95% CI = 1.05–20.85, P = 0.042) (Fig. 4A-4D). What's more, metastatic or recurrent patients had higher risk scores than those without metastasis or recurrence. (Fig. 4E and 4F).

Fig. 4

Relationships between apoptosis-associated gene signature and clinical parameters in GSE21257 and GSE39055. (A and B). Univariate and multivariate cox regression analyses revealed the risk score was the independent prognostic factor in GSE21257. (C and D). Univariate and multivariate cox regression analyses validated risk score and recurrence were two independent prognostic factors in GSE39055. (E). The correlation between apoptosis-associated risk score and osteosarcoma metastasis in GSE21257. (F). The relationship between apoptosis-associated risk score and osteosarcoma recurrence in GSE39055.

Then the correlations between clinical parameters and individual genes in signature were explored. In GSE21257, the expressions of HSPB1 and IER3 were not associated with the age of diagnosis or gender (Supplementary Figure 2A,2B,2E,2F). On the contrary, these two genes were significantly related to metastasis and apoptosis-associated risk (Supplementary Figure 2C,2D,2 G,2H). Similar results were obtained in GSE39055 except that HSPB1 was not correlated with patients’ recurrence status (Supplementary Figure 3A–3H).

Construction of the nomogram

In order to accurately estimate the prognosis of osteosarcoma patients, the nomogram was established based on the clinical variables and risk scores. The relationships between nomogram scores and patients’ 3- or 5-year survival rates were shown in Fig. 5A and Supplementary Figure 4A. The 3-year and 5-year calibration curves revealed the high predictability of nomograms and the C-index values (0.86 and 0.95) indicated the nomogram was reliable (Fig. 5B, 5C and Supplementary Figure 4B,4C).

Fig. 5

Establishment of the nomogram in the GSE21257. (A) The nomogram to predict the 3- and 5-year survival risk of osteosarcoma patients. (B) The calibration curve of the 3-year survival. (C) The calibration curve of the 5-year survival.

Gene set enrichment analysis of the apoptosis-associated prognostic signature

To identify the dysregulated pathways related to the apoptotic risk, GSEA was performed. The results showed that in GSE21257 and GSE39055 (Supplementary Figure 5), high-risk patients were accompanied by numerous downregulated immune functions such as natural killer cell mediated cytotoxicity, B cell receptor signaling pathway and intestinal immune network for IgA production. These results indicated that there might be potential close relationships between apoptosis and immune system.

Exploration of the correlation between apoptosis and immune infiltration

The immune infiltration status of 22 immune cells in osteosarcoma was evaluated through CIBERSORT algorithm. The abundance ratios of osteosarcoma immune cells with statistically significant were shown in Fig. 6A and Supplementary Figure 6A. The interrelations between these cells were displayed in Fig. 6B and Supplementary Figure 6B Then the specific relations between apoptosis-associated risk and immune infiltration were conducted in ssGSEA analysis. Interestingly, almost all immune cells (14) and functions (13) were dysregulated in GSE21257 (Fig. 6C, 6D). Meanwhile, two immune cells and four immune functions were differentially expressed in GSE39055 (Supplementary Figure 6C, 6D). The overlapped items included one immune cell (Th1 cells) and four immune functions (CCR, Check-point, MHC class I, and T cell co-inhibition), while MHC class I showed the opposite trend and was removed (Fig. 6E).

Fig. 6

The relationships between apoptosis and immune infiltration in GSE21257. (A) The infiltrative proportion of 22 immune cells in osteosarcoma with statistically significant. (B) The heatmap of correlations between different immune cells in osteosarcoma. The number and circles in each tiny box represented the corresponding correlation value between two cells. The blue color indicated positively related and the red indicated negatively related. (C) The correlations between risk score and different immune cells. (D) The association between risk score and different immune features. (E) The overlapped immune cells and functions among GSE21257 and GSE39055. The red color implied the opposite trend in two datasets. NS, not significant; *P < 0.05; **P < 0.01; ***P < 0.001.

Validation of the apoptosis-associated genes in osteosarcoma samples

We screened out eight patients who were diagnosed as osteosarcoma by pathology in Ruijin Hospital and separated them into non-metastasis (4) and metastasis (4) groups. Then the expressions of apoptosis-associated genes in human osteosarcoma samples were validated by IHC assays. The results revealed that HSPB1 was down-regulated in non-metastasis patients (Supplementary Figure 7A) while IER3 showed the opposite trends (Supplementary Figure 7B).

Discussion

Osteosarcoma was one of the most common primary bone malignant tumor which often developed metastasis [19]. Nowadays, most visible osteosarcoma could be resected by surgery and some micro-satellite lesions were able to be eliminated by chemotherapy or targeted therapy. However, many patients still suffered from metastasis after radical treatment and their prognoses were extremely poor [20]. The mechanism of cancer metastasis was first proposed by Paget in 1889 which was widely known as ‘seed and soil’ hypothesis [21]. After that, various mechanisms related to tumor metastasis were revealed, such as enhancement of angiogenesis, epithelial-mesenchymal transition (EMT), chemoresistance and so on [22]. Among these hallmarks, apoptosis was intensively studied and regarded as one of the most important prognostic markers.

Apoptosis, one of the widely known patterns of cell death, was related to numerous physiological and pathological alterations in cells. The whole process of apoptosis was triggered by two pathways. One is the extrinsic pathway which was mainly activated by death receptors at the cell surface. The other one was the intrinsic pathway which was mainly induced by drugs, endoplasmic reticulum (ER) stress, perforin and so on [23]. These pathways could activate effector caspases to eliminate the inhibitor of apoptosis. Then the activation of caspase-activated deoxyribonuclease (CAD) would alter the structure of DNA and dysregulate the activity of cell structure-associated proteins. After that, the cell integrity was destructed and the cells were broken into apoptotic bodies [24]. During the occurrence and development of malignant tumors, the imbalance between proliferation and apoptosis always existed. In the last few decades, numerous apoptosis-associated genes had been identified and some of them were related to tumor progression in various cancers [25]. In osteosarcoma, these genes also played important roles. For instance, Lu et al. found that the transcriptional repressor gene ZBTB7A could inhibit the expression of lncRNA GAS5 and suppress the ER stress-induced cell apoptosis, which promoted the progression of osteosarcoma [26]. In addition, some researchers also found that these apoptosis-associated genes could be combined as the signature to predict the prognosis of patients [27,28]. However, the correlation between osteosarcoma and apoptosis-associated prognostic signature was unclear.

In this study, 464 DEmRNAs associated with osteosarcoma metastasis were identified in GSE21257. Then nine genes were screened out through the intersection between these DEmRNAs and 161 apoptosis-associated genes in MSigDB. The apoptosis-associated prognostic signature was constructed through univariate and multivariate cox regression analyses. After that, all patients in GSE21257 were separated into high- and low-risk groups based on the median value of risk score. The KM survival curve showed that the apoptosis-associated risks were negatively related to the survival rates of patients and the AUC values of ROC curves were over 0.7, which indicated the robustness of this prognostic signature. What's more, the risk score was regarded as the independent prognostic factor and its expression was increased in patients with metastasis.

The nomogram was a powerful tool that helped clinicians to estimate the prognosis of patients. In addition, patients with high nomogram evaluative scores could receive more aggressive treatment in the early stage. In our research, the clinical variables and risk scores were integrated to construct the nomogram. The C-index and calibration curves revealed the reliability of this apoptosis-associated prognostic signature. The same results were obtained in another dataset GSE39055. To sum up, we thought that the apoptosis-associated signature was a convincible predictor for osteosarcoma prognosis and could be applied in clinical.

The tumor microenvironment mainly comprised various types of cells (endothelial, fibroblastic, immune cells, etc.) and extracellular components [29]. In the process of tumor apoptosis, the components of tumor microenvironment such as immune cells were often dysregulated. For example, Young et al. revealed that the enhancement of autophagy could inhibit TNFɑ-induced apoptosis and thus protected tumors from T cell-mediated cytotoxicity [30]. Sun et al. found that in non-small cell lung cancer, elevated LINC00301 suppressed cell apoptosis and promoted the infiltration of tumor-suppressive immune cells [31]. Numerous studies also reported the same results in osteosarcoma. For instance, Zou et al. found that C3AR1 promoted the apoptosis of osteosarcoma cells and was closely related to the infiltration of several immune cells such as macrophages M1, macrophages M2, CD8+ T cells and so on [32]. In our research, the GSEA results revealed that the apoptosis-associated risk related to several immune pathways such as B cell receptor signaling pathway and intestinal immune network for IgA production. Next, the CIBERSORT algorithm and ssGSEA were conducted to evaluate the immune infiltration of osteosarcoma. The results indicated that high apoptosis-associated risk scores correlated with decreased infiltrative immune cells (Th1 cells) and inhibitory immune functions (CCR, Check-point, T cell co-inhibition). These results were also validated in other studies. For example, in ovarian cancer, Chen et al. revealed that artesunate could enhance the apoptosis of cancer cells by increasing the proportion of Th1 cells [33]. In a word, our research revealed that the apoptosis might influence the immune infiltration of osteosarcoma and thus promoted the progression and metastasis of cancer.

What's more, the expression of individual mRNAs in the signature was investigated. For IHC assays, the HSPB1 was up-regulated in osteosarcoma patients with metastasis while IER3 exhibited the contrary tendency. Numerous studies also revealed that the mRNAs in our apoptosis-associated signature played important roles in tumor development. HSPB1, also known as HSP27, was first purified and studied from 1987 to 1988 [34,35]. It was a member of heat-shock protein family and was expressed when cells experienced external or internal stresses such as heat shock, oxidative stress, hypertonic stress and so on [36]. It was a repressor of apoptosis by inhibiting both Fas ligand pathway and mitochondrial pathway. It could also active molecular chaperones, inhibit cytoskeleton destruction and regulate intracellular redox action, which eventually protected cells from stresses [37]. The correlation between osteosarcoma and HSPB1 was first mentioned by Hiroshi in 1996. They found that overexpression of HSPB1 was related to poor preoperative chemotherapy effects and it could be regarded as the independent prognostic factor [38]. After that, various studies also revealed that HSPB1 was related to distant metastasis, chemoresistance, autophagy and poor prognosis in osteosarcoma [39], [40], [41]. IER3, also known as IEX-1, was an immediate early response gene which was mainly induced by cellular disturbance such as DNA damage, growth factors, cytokines and viral infection [42]. It was widely expressed in human tissues and was involved in various cellular processes such as apoptosis, proliferation, differentiation and so on [43], [44], [45]. Its relationship with tumor initiation and development was widely investigated. Some studies revealed that overexpression of IER3 would enhance tumor cells’ sensitivity to radiotherapy as well as chemotherapy and promoted the apoptosis of cancers [46,47]. Zhang et al. found that overexpression of IER3 would decrease the apoptosis of activated T cells and increase the effective state of immune response [48]. However, their specific role in osteosarcoma was unclear and should be investigated in the future.

To our knowledge, this was the first study that investigated the relationship between apoptosis-associated gene signature and osteosarcoma patients’ prognosis. However, our research still had some limitations. First, the sample size of dataset was relatively small. This might relate to the low incidence of osteosarcoma and lack of large-scale studies. Second, in consideration of the complicated processes of apoptosis, some genes related to apoptosis might be omitted. Last but not least, our study was a retrospective research and this signature should be validated in other prospective studies or clinical trials in the future.

Conclusions

In conclusion, we identified and validated the apoptosis-associated prognostic signature related to osteosarcoma metastasis. What's more, we revealed that this signature was related to the immune infiltration of osteosarcoma. This gene signature could act as the biomarker to predict the prognosis of patients and might provide novel therapeutic targets for clinical therapy. Further studies and clinical trials should be conducted to clarify the specific mechanisms of our findings.

CRediT authorship contribution statement

Yucheng Fu: Conceptualization, Methodology. Zhijian Jin: Conceptualization, Methodology. Yuhui Shen: . Zhusheng Zhang: Data curation, Writing – original draft. Meng Li: Software, Validation. Zhuochao Liu: Software, Validation. Guoyu He: Data curation, Writing – original draft. Jintao Wu: Data curation, Writing – original draft. Junxiang Wen: Visualization, Investigation. Qiyuan Bao: Visualization, Investigation. Jun Wang: Writing – review & editing. Weibin Zhang: Supervision, Writing – review & editing.

Conflict of interests

None.