Introduction

Osteosarcoma (OS) is the most common primary malignant tumor in children and adolescents. However, its incidence rate is quite low, with approximately 3610 new cases diagnosed annually worldwide1,2. Although the mortality rate of OS has dramatically decreased due to neoadjuvant multiagent systemic chemotherapy, advanced surgical techniques, and precision radiotherapy, the survival rate is still not satisfactory, especially for those patients with tumor metastasis and recurrence3,4. It has been challenging to improve the prognosis of OS.

Numerous scientists have studied the gene mechanisms and treatments of OS. The development of humanized in vitro and humanized mouse models with similar tumor microenvironments has greatly facilitated related research5,6. In vitro 3D tumor models showed that MAPK, TGFβ/SMAD, PI3K/AKT, JAK/STAT, Notch and Hedgehog signaling transition molecules exhibit significantly increased expression7. Nigris et al. also reported that OS cells cultured in scaffolds showed a dramatic increase in angiogenic factors such as PDGFB, TGFB1, VEGFA and VEGFB8. HER2, CXCR4 and HIF2α were found to be typical biomarkers involved in tumor growth and the homing of cancer cells to distant sites9. Humanized in vitro models of osteosarcoma also showed improved drug resistance to doxorubicin compared with scaffold-free spheroids10. Regarding treatment, Dimas et al. found that targeting the farnesylation of Ras by self-assembling nanoparticles encapsulating zoledronic acid simultaneously inhibited the Ras/ERK1/2/HIF-1α and Ras/Akt/mTOR axes, leading to reduced growth and increased intratumour necro-apoptosis and T-cell infiltration of chemo-immune-resistant osteosarcomas11. Targeting mitochondrial complex I prevented HIF1-MIF activation, leading to TAM accumulation and vascular architecture remodeling and ultimately resulting in the inhibition of osteosarcomas grown in the humanized model12. Immunotherapy drugs such as nivolumab have also been investigated in a humanized mouse model by Zheng et al., and the results showed that the nivolumab-treated group exhibited a similar primary tumor volume and growth rate but a significantly lower lung metastasis rate than the control group. Further analysis showed that CD4 + and CD8 + lymphocytes were more frequently observed in the lungs of the nivolumab-treated group than in the lungs of the control group, but no statistically significant differences were observed in the primary tumors from both groups13. All these results indicated that immunotherapy targeting the tumor microenvironment may be a promising strategy for osteosarcoma treatment.

The immune microenvironment of OS has been investigated since the accidental discovery of tumor responses following bacterial toxin treatment. Its functions are diverse and complex and have not been fully understood until now14,15. Previous results showed that osteosarcomas are infiltrated mainly by macrophages and T cells16,17,18. Studies have shown that high infiltration levels of macrophages and CD8 + T cells are associated with reduced metastasis and improved survival in OS19,20,21. In contrast, high infiltration levels of antigen-presenting cells have been reported to lead to unfavorable outcomes22. In recent years, several studies have focused on the detailed molecular mechanisms; for example, Wang et al. found that OS tumor cells could release PD-L1 to promote the metastatic process by inhibiting the immune response23. Troyer et al. reported that indoleamine dioxygenase could inhibit dendritic cells from producing neoantigens, thus indirectly leading to immune escape24. A previous study also reported that all-trans retinoic acid could inhibit M2 polarization of macrophages to repress the OS lung metastatic process25. In addition, the VEGF, IL-10A, TGF-β, and STAT3 pathways have also been found to facilitate the immunosuppressive microenvironment by influencing bone marrow-derived suppressor cells, macrophages and stromal fibroblasts15. Together, these data highlight the important role of the immune microenvironment in patients with OS. However, the mechanism by which infiltrating immune cells regulate the pathogenesis and development of OS remains largely unknown. In addition, it is not ideal to predict the prognosis and therapeutic response solely using infiltrating immune cells of OS patients. With the appearance of high-dimensional datasets and advanced bioinformatics algorithms26,27,28, it is now possible to analyze the comprehensive interactions between biological phenotypes and the tumor immune microenvironment, thus facilitating research on the molecular characteristics affecting immune cell infiltration, the response to immunotherapy, and the prognosis of OS patients.

In our study, we gained a deeper understanding of the OS immune microenvironment by utilizing OS cohorts from four public databases, which showed that M0 macrophages were the most consistently infiltrating immune cell type. We then screened for genes correlated with M0 macrophages in the Therapeutically Applicable Research to Generate Effective Treatments (TARGET) and GSE21257 datasets under the supposition that dominant immune cells have more vital effects on OS prognosis. A prognostic immune signature was built by the least absolute shrinkage and selection operator (LASSO) Cox model using the abovementioned M0 macrophage-correlated genes from the TARGET training cohort. The predictive values of this model were further validated in two other independent testing cohorts and verified by comparing them to previous prognostic models. Then, we attempted to identify the related signaling pathways, intrinsic molecular subtypes, hot immune-targeted gene expression, and distributions of immune cells between each risk score subtype. Furthermore, the relationship between the risk score and predicted chemotherapy sensitivity was also assessed.

Methods

Data collection

Gene chip expression data and the clinical information of OS patients from the TARGET database were obtained from the website https://ocg.cancer.gov/programs/target/data-matrix on May 1, 2022. Osteosarcoma was used as the only key word to select suitable datasets published up to May 2022 for our study on the Gene Expression Omnibus (GEO) datasets website (https://www.ncbi.nlm.nih.gov/gds/). Then, the datasets were searched by the following inclusion criteria: (1) the diagnoses were pathologically confirmed; (2) the dataset had Homo sapiens samples; (3) the platform contained whole-genome information; (4) the dataset included patient survival data; and (5) the sample size of the datasets was more than 30. After that, only the GSE21257, GSE16091, and GSE39055 datasets were included in our study. Therefore, the normalized mRNA expression data and survival data of the above datasets were downloaded from GEO (https://www.ncbi.nlm.nih.gov/geo/query/acc.cgi?acc=GSE21257, https://www.ncbi.nlm.nih.gov/geo/query/acc.cgi?acc=GSE16091, https://www.ncbi.nlm.nih.gov/geo/query/acc.cgi?acc=GSE39055).

If the same patient provided two or more tumor samples to those datasets, only the data corresponding to the primary lesion were selected according to the sample numbers. All patients were included in our study. Gene symbols were extracted from the provided documents from those dataset websites. We conducted signal intensity normalization across the arrays of the above datasets using the normalizeBetweenArrays function from the limma package in R software.

Tumor immune cell infiltration calculation

The infiltration levels of 22 immune cells in the TARGET, GSE21257, GSE16091, and GSE39055 datasets were quantified by the “CIBERSORT” R package using the 1,000 permutations method. The reference immune cell signatures were downloaded from the Supplementary Information on the article’s website (https://www.nature.com/articles/nmeth.3337#MOESM207)29,30,31. The ESTIMATE score, stromal score, and immune score of each patient were assessed by the “ESTIMATE” R package using the gene expression data from each patient sample32. All parameters in the R equation throughout the calculation process were set to default.

M0 macrophage correlated gene screening

This section attempts to select genes associated with M0 macrophages. Pearson's correlation analysis was used to search for M0 macrophage-correlated genes with cutoffs of |correlation coefficient|> 0.3 and P value < 0.05 in the TARGET and GSE21257 datasets, which have a relatively large sample size. Next, we merged those correlated genes from the two datasets to obtain a full correlated gene list.

Gene signature construction and validation

Gene signatures were generated by inputting M0 macrophage-correlated genes into the LASSO Cox regression model using the TARGET dataset. The “glmnet” package in R was used to complete the regression process. After the genes were selected, a multivariate Cox analysis was used to calculate the corresponding coefficients. The signature risk score was calculated by the sum of the products of each gene and its corresponding coefficients as follows: score = (CPE*0.7021) + (FHL5*− 0.9794) + (GBP1*− 0.4718) + (GNLY*0.9393) + (GPR82*− 1.8384) + (IL18RAP*− 1.9630) + (LILRA2*0.4310) + (NDRG4*− 0.7661) + (PLB1*− 0.6822). Moreover, the risk scores for the GSE21257 and GSE39055 datasets were also calculated using the same method with the coefficients derived from the TCGA dataset for validation. OS patients from each dataset were divided into two groups (low- and high-risk groups) according to the median value of the risk score to maximize the statistical power and provide the most reliable results based on a fixed sample size in those databases33. An overall survival analysis between the two groups was conducted using Kaplan–Meier curves with Wilcoxon log-rank tests. Multivariate Cox regression models were applied to test the independent predictive value of our signature. The prognostic accuracy of the risk score in different datasets was determined using the time-dependent receiver operating characteristic curve (ROC) and incident/dynamic (I/D) area under the curve (AUC) through the “timeROC” and “risksetROC” packages in R separately34.

Gene signature comparison

We screened the studies using the keyword “gene signature prognosis osteosarcoma” on PubMed. In addition, the inclusion criteria were as follows: (1) the journal impact factor was more than 5 and had a good reputation; (2) the online publication date range was from March 01, 2020, to March 01, 2022; and (3) the gene signature was constructed by messenger RNA to be easily validated in other datasets. After the studies were selected, the gene signatures were recalculated by a multivariate Cox proportional hazards model using data in a training cohort. Kaplan–Meier curves, time-dependent ROCs and I/D AUCs of the gene signatures using validation cohorts were chosen to compare the prognostic accuracy.

Gene set enrichment analysis (GSEA)

GSEA was launched to study the biological functions of our gene signature in the TARGET, GSE21257, and GSE39055 datasets using the “clusterProfiler” R package based on the GO and KEGG analyses individually. The GO functions and KEGG pathways with P values less than 0.05 in each dataset were merged. Finally, the biological functions with consistent positive or negative values of enrichmentScore were reserved to find the real difference of GO functions and KEGG pathways between low- and high-risk groups discriminated by our signatures.

Correlation between signature risk score and chemotherapy sensitivity

The IC50 values of chemotherapy drugs (bleomycin, cisplatin, doxorubicin, etoposide, and methotrexate) in each osteosarcoma sample in the TARGET, GSE21257, and GSE39055 datasets were accessed through the “pRRophetic” package in R35, which was built based on the Genomics of Drug Sensitivity in Cancer database (www.cancerRxgene.org)36. Then, the samples were classified into low- or high-risk groups by the median of our gene signature risk score. Finally, the IC50 values of chemotherapy drugs between different groups were analyzed by the Mann‒Whitney rank test. A P value < 0.05 was defined as statistically significant.

Statistical analysis

We conducted all statistical analyses with R version 4.0.5 (R Foundation for Statistical Computing, Vienna, Austria) and GraphPad Prism 6.01 (GraphPad Software, Inc., San Diego, CA, USA). OS patients were divided into two groups by the median of our signature risk score in the TARGET, GSE21257, and GSE39055 datasets. Violin plots were performed in Hiplot (https://hiplot.com.cn), a comprehensive web platform for scientific data visualization. The comparison of 22 immune cells, immunotherapy-targeted genes, ESTIMATE scores, stromal scores, and immune scores between low- and high-risk groups were analyzed using a Mann‒Whitney rank test. Unless otherwise specified, a P value < 0.05 was defined as statistically significant.

Ethics approval

Our study is based on open-source data (TRGET and GEO). Ethical review and approval were not required for the study on human participants in accordance with the local legislation and institutional requirements.

Consent to participate

All methods were carried out in accordance with the relevant guidelines and regulations.

Results

Clinical characteristics of the studied datasets

Figure 1 shows the flowchart of this study. The clinical characteristics of OS patients in the TARGET, GSE21257, GSE16091, and GSE39055 datasets are shown in Table 1. In the TARGET and GSE21257 datasets, the age, sex, tumor location, and metastatic status at diagnosis were similar. The GSE21257 dataset contained the most detailed clinical information except for radiotherapy treatment and had a similar Huvos grade distribution as the GSE39055 dataset. GSEGSE16091 only included survival information. The median overall time and survival curves (Supplemental Figure 1) were close to each other in all datasets.

Figure 1
figure 1

A schematic diagram of dataset analysis.

Full size image
Table 1 Characteristics of OS patients in the TARGET, GSE21257, GSE16091 and GSE39055 datasets.
Full size table

Estimation of tumor immune cell infiltration in osteosarcoma

To obtain the landscape of tumor immune cell infiltration, we applied the CIBERSORT algorithm to the selected datasets. As shown in Fig. 2, the major immune cells in the TARGET and GSE16091 datasets were CD8 + T cells and M0 and M2 macrophages, while the major immune cells in the GSE21257 dataset were M0 and M2 macrophages. For the GSE39055 dataset, M0 macrophages were the most dominant infiltrating immune cells. Other immune cell infiltration levels were quite low in all datasets. In general, M0 macrophages were the dominant immune cells across all datasets. As shown in Supplemental Figure 1, the survival difference between each dataset was not significant, and we believe that M0 macrophages are influential on the overall survival of osteosarcoma. Therefore, the M0 macrophage-associated genes were screened by the methods described in the “Methods and Materials”, and the selected genes are shown in Supplemental Table 1.

Figure 2
figure 2

Estimation of tumor immune cell infiltration in osteosarcoma by the CIBERSORT algorithm. The percentage of infiltration of 22 immune cell types in the TARGET (A), GSE21257 (B), GSE16091 (C), and GSE39055 (D) datasets showed that M0 macrophages were the dominant immune cell type.

Full size image

Construction and validation of prognostic gene signature for osteosarcoma

The above selected genes were used to build a prognostic signature in OS patients through LASSO Cox regression analysis by the TARGET dataset with 86 patients as the discovery cohort (Fig. 3A). An optimal 9-gene prognostic signature was made (Fig. 3B,C)). The biological functions and the coefficients of signature genes are shown in Table 2. The signature risk scores are equal to the sum of the product of the expression value and coefficient of each gene. We chose only the GSE21257 and GSE39055 datasets as the validation cohort to measure the prognostic value of the signature-based risk score, and the GSE16091 dataset was excluded owing to the absence of gene expression data for GPR82 and PLB1. We used the log-rank test to study the association of the risk score with survival in the TARGET dataset. As expected, high-risk patients had significantly shorter survival than low-risk patients (log-rank test, P < 0.0001) (Fig. 3D). The same tendency was confirmed in the GSE21257 and GSE39055 datasets (log-rank test, P = 0.034 and P = 0.031) (Fig. 4A,B). In addition, the AUCs of the time-dependent ROCs including the signature risk score and selected clinical factors (age, sex, and Huvos grade) reached 0.793, 0.885, and 0.839 for 1-, 3-, and 5-year OS in the TARGET dataset, respectively, which were higher than the AUCs of the time-dependent ROCs including only clinical factors (Fig. 5A, B, C). Similar results were also observed in the GSE21257 and GSE39055 datasets (Fig. 5E,F,G,I,J,K).

Figure 3
figure 3

Construction of a prognostic gene signature. (A) Venn diagram of the M0 macrophage-associated genes in the TARGET and GSE21257 datasets. The numbers in each area represent the gene numbers in each group. (B) Cross-validation for tuning parameter screening upon LASSO regression analysis. (C) Screening of the optimal parameter (lambda) at which the vertical lines were drawn. (D) Kaplan‒Meier overall survival analysis of the gene signature risk score in OS of the TARGET dataset. (E) Distribution of signature gene expression profiles along with survival status in different signature risk score groups in TARGET datasets. (F) Forest plot showing the results of multiple factors in the Cox regression analysis of the gene signature risk score with other clinical characteristics in OS of the TARGET dataset.

Full size image
Table 2 Prognostic genes obtained from the LASSO Cox regression model.
Full size table
Figure 4
figure 4

Validation of the nine-gene prognostic signature. Kaplan‒Meier overall survival analysis of the gene signature risk score in OS from the GSE21257 (A) and GSE39055 (B) datasets. Distribution of signature gene expression profiles along with survival status in different signature risk score groups in the GSE21257 (C) and GSE39055 (D) datasets. Forest plot shows the results of multiple factors in the Cox regression analysis of the gene signature risk score with other clinical characteristics in overall survival from the GSE21257 (E) and GSE39055 (F) datasets.

Full size image
Figure 5
figure 5

Prognostic values of the signature risk score in the training and validation cohorts. Time-dependent ROCs at 1 (A), 3 (B), and 5 (C) years and I/D AUC (D) for OS in the TARGET dataset. Time-dependent ROCs at 1 (E), 3 (F), and 5 (G) years and I/D AUC (H) for OS in the GSE21257 dataset. Time-dependent ROCs at 1 (I), 3 (J), and 5 (K) years and I/D AUC (L) for OS in the GSE39055 dataset. ROC, receiver operating characteristic curve; I/D AUC, incident/dynamic area under the curve.

Full size image

We compared the C index of I/D AUC to assess the discriminative accuracy of the prediction model with and without the signature risk score. The results showed that the C index increased to 0.831 (95% CI, 0.750–0.868 P < 0.0001), 0.690 (95% CI, 0.527–0.754; P < 0.0001), and 0.775 (95% CI, 0.584–0.840; P < 0.0001) when the signature risk score was added to the prognostic model in the TARGET, GSE21257, and GSE39055 datasets (Fig. 5D,H,L). Figure 3E and Fig. 4C,D show the distribution of signature gene expression profiles and survival statuses in different risk score groups for OS in all datasets. Figure 3F and Fig. 4E,F show that the high-risk score was significantly associated with a risk of increased mortality in OS patients in the multivariate Cox regression models after adjusting for age, sex, and Huvos grade in all datasets (HR = 51.871, 3.204 and 6.663, 95% CI = 6.898–391.183, 1.066–9.637 and 1.227–36.171, P < 0.001, = 0.038 and = 0.028, respectively), which indicates that the signature risk score is an independent prognostic factor for OS patients.

Investigation into the differences in immune cells, immunotherapy-targeted genes, and immune scores between high- and low-risk groups

As the signature risk score reflects the tumor immune activity of OS, we further investigated the differences in immune cells, immunotherapy-targeted genes, and immune scores between different risk groups from all datasets. The immune cell infiltration results were inconsistent between the three datasets (Supplemental Figure 2). M0 macrophages were found to be higher in high-risk groups than in low-risk groups in the TARGET and GSE21257 datasets (Fig. 6A, B). Activated dendritic cells were found to be lower in the high-risk groups than in the low-risk groups in the GSE21257 and GSE35099 datasets (Fig. 6B, C). For the results of six hot immunotherapy-targeted genes and immune scores, no significant differences were found in the GSE35099 datasets (6F, I), which may be due to its relatively small sample size. TIM-3, LAG-3 and all immune scores were all found to be lower in the high-risk groups than in the low-risk groups in the TARGET and GSE21257 datasets (Fig. 6D,E,G,H). In summary, our signature risk score may reflect the immune activity of osteosarcoma. In addition, LAG3 may be the most effective immunotherapy target because its expression was relatively high in all datasets (Fig. 6D,E,F). Owing to the low incidence rate of OS, the OS datasets did not contain large sample sizes. Therefore, adjustments for multiple testing were not made to reserve important clues.

Figure 6
figure 6

Association between immune cell infiltration levels, immunotherapy-targeted gene expression, ESTIMATE immune scores, and signature risk score of OS. The immune cell infiltration types found to be significantly associated with the signature risk score of OS in the TARGET (A), GSE21257 (B), and GSE39055 (C) datasets. Association between immunotherapy-targeted gene expression and the signature risk score of OS in the TARGET (D), GSE21257 (E), and GSE39055 (F) datasets. Association between ESTIMATE immune scores and signature risk scores of OS in the TARGET (G), GSE21257 (H), and GSE39055 (I) datasets.

Full size image

Comparison of gene signatures

Thirteen previously reported studies were included in the comparison after the selection process, which are shown in Table 3. Kaplan–Meier curves (Fig. 7) and univariate Cox models (Table 4) showed that only the difference in overall survival for the nine-gene signature model was significant in the training and testing cohorts. Additionally, the AUCs of the time-dependent ROCs of the signature risk score were also found to be greater than the AUCs of the time-dependent ROCs of previous signatures in the TARGET and GSE39055 datasets (Fig. 8). The results showed that our nine-gene prognostic signature is more robust than previously reported gene signatures.

Table 3 Candidate research for comparison to our signature.
Full size table
Figure 7
figure 7

Comparisons of the gene signature with previously published gene signatures in the TARGET, GSE21257, and GSE39055 datasets using the Kaplan–Meier estimator. The results are in bold and considered significant if the P value < 0.05.

Full size image
Table 4 Comparison of the nine-gene signature with previous published models in a univariate Cox analysis.
Full size table
Figure 8
figure 8

Comparisons of the gene signature with previously published gene signatures in the TARGET, GSE21257, and GSE39055 datasets using time-dependent ROC curves.

Full size image

Identification of signature-associated biological functions and pathways in OS

The tumor sample tissues of the three datasets were dichotomized into high- and low-risk groups according to the median of the nine-gene signature risk score. A GSEA was performed to identify the signature-associated biological functions and signaling pathways, and then the results from all datasets were merged to attain reliable results. We found that the olfactory transduction receptor signaling pathway and its related biological functions were downregulated in high-risk groups compared to low-risk groups (Fig. 9). Antigen processing- and presentation-related biological functions were also found to be upregulated in high-risk groups compared with low-risk groups. However, none of the related signaling pathways were found to be significantly correlated with different risk groups. These results indicate that the olfactory transduction receptor signaling pathway may play a role in the difference between signature-predicted outcomes.

Figure 9
figure 9

Gene set enrichment analysis of GO and KEGG pathways in OS between different signature risk score groups in the training and validation cohorts. The results of GO functions between different signature risk score groups in the TARGET (A), GSE21257 (B), and GSE39055 (C) datasets. The results of KEGG pathway analysis between different signature risk score groups in the TARGET (D), GSE21257 (E), and GSE39055 (F) datasets.

Full size image

Signature risk scores and chemotherapy sensitivity

To investigate whether the immune-correlated gene signature is an independent prognostic factor of OS patients, we need to adjust it based on previously reported important prognostic factors, including histologic response to chemotherapy and presence of metastases. Supplemental Figure 3 shows the results and indicates that our signature is an independent OS prognostic factor. However, not all target patients received neoadjuvant chemotherapy; thus, this group of patients lacked Huvos grade information. To validate the independence of our signature, we explored the relationship between the risk score and chemotherapy sensitivity. IC50 was calculated using the “pRRophetic” R package to predict the treatment response to chemotherapy drugs (bleomycin, cisplatin, doxorubicin, etoposide, and methotrexate). As we anticipated, there was no significant difference in the sensitivity of all drugs in the low- and high-risk groups (Fig. 10). The results are consistent with the multivariate Cox analysis. GSEA also did not reveal any biological functions or pathways responsible for chemotherapy resistance. Therefore, the overall survival difference predicted by our signature is probably due to the various backgrounds of the immune microenvironment of OS.

Figure 10
figure 10

Gene signature risk score in the role of chemotherapy. (A) The correlation between different signature risk scores of OS patients and estimated IC50 values of bleomycin, cisplatin, doxorubicin, etoposide, and methotrexate in the TARGET (A), GSE21257 (B), and GSE39055 (C) datasets. The Wilcoxon signed-rank test was used to compare the estimated IC50 value of chemotherapy between groups.

Full size image

Discussion

Although OS is the most common primary bone cancer in children and young adults, it is a very rare cancer, with approximately 400 new cases diagnosed annually in the USA37. The incidence peaks in adolescence and in old age38. The most common early symptom is ostalgia, which is easily confused with growing pain. The most common location of OS is in the metaphysis around the knee joint, followed by the proximal tibia and humerus. The majority of patients are therefore diagnosed with localized disease. Among the patients in the progression stage, the lungs and other bones are the most common metastatic locations39. The prognostic factors include tumor site, tumor size, tumor resectability, histologic response to chemotherapy, and presence of metastases40,41,42. Chemotherapy consisting of high-dose methotrexate and other drugs before and after definitive surgery has been used as the standard treatment strategy for localized osteosarcoma over the past decades43,44,45,46. In patients with relapsed and/or metastatic osteosarcoma, surgery and chemotherapy have some effect3,47,48,49. However, the prognosis of those patients was a median event-free survival of less than 4 months, which is unsatisfactory50. An increased understanding of the disease from well-annotated tissue banks through newly developed technologies has revealed its heterogeneity and molecular aberrations, which offer new insight into targeted therapy and immunotherapy.

Recently, immunotherapy has facilitated a revolution in many kinds of solid malignant tumor treatments. Although its efficiency is not so high, it can lead to considerable survival benefits among responders51. Osteosarcoma humanized in vitro and humanized mouse models with similar tumor microenvironments in several previous studies have indicated that immunotherapy targeting the tumor microenvironment may be a promising strategy for osteosarcoma treatment11,12,13. Studies aiming at finding the precise immune response population have also suggested that the tumor immune microenvironment is the critical factor52. Therefore, it is essential to better understand and characterize the immune status of OS for the precise prediction of the immune response and the development of immunotherapy. The OS immune microenvironment consists of a network of immune cells that function as ideal grounds for tumor proliferation and progression15,53.

Our study first found that M0 macrophages were the dominant immune cells in the TARGET and GEO databases using the CIBERSORT algorithm. However, the direct association between M0 macrophage infiltration and the overall survival of OS patients was not consistent in different independent datasets (Supplemental Figure 4). Based on the hypothesis that the dominant immune cells play the most important role in the OS microenvironment, we screened for genes associated with M0 macrophages and built a new nine-gene signature for OS prognosis using clinical information from the TARGET dataset. In addition, the prognostic signature was proven to be effective and reliable in validation cohorts consisting of two independent GEO datasets. Previously published prognostic gene signatures were compared with the newly discovered gene signature using publicly available datasets. The results showed that no signatures besides ours offer a precise prediction value of OS survival in all datasets, indicating the priority of the nine-gene signature. Multivariate Cox regression models indicate that the signature risk score is an independent prognostic factor for OS patients regardless of age, sex, or Huvos grade in all datasets. Furthermore, the association between the signature risk score and chemotherapy sensitivity also showed that there was no significant difference in the sensitivity of all drugs between the low- and high-risk groups. Therefore, the different outcomes predicted by our signature are probably due to the various immune microenvironments of OS. To our knowledge, this is the first time a study has been conducted on the dominant infiltrating immune cell profiles in several datasets to find similarities between OS tumor immune microenvironments. As we expected, a robust gene signature was generated based on the discovered similarities.

M0 macrophages are slightly elongated unstimulated macrophages, which are considered to be theoretically inactivated54. Previous studies have suggested that M0 macrophages are associated with poor prognosis and metastatic disease in lung adenocarcinoma55,56, hepatocellular carcinoma57, osteosarcoma53, pancreatic adenocarcinoma58, melanoma59, colorectal cancer60, gastric cancer61, and glioblastoma62. Recent studies have shown that new tumor-infiltrated M0 macrophages induce pancreatic cancer cell death via TNF-α secretion, while M1, M2, and tumor-associated macrophages do not harbor antitumorigenic activities63. However, mesenchymal stem cells could activate STAT6 and induce M2 polarization through disease progression, leading to anti-inflammatory functions64, which may explain the relationship between high M0 macrophage infiltration and poor prognosis.

Our prognostic signature consisted of nine genes (Table 2): CPE, FHL5, GBP1, GNLY, GPR82, IL18RAP, LILRA2, NDRG4, and PLB1. In our research, CPE was determined to be an unfavorable biomarker for OS prognosis, while PLB1 exhibited protective effects in OS (Supplemental Figure 5). CPE encodes carboxypeptidase E, which belongs to the carboxypeptidase family and is reported to be involved in the biosynthesis of hormone and neurotransmitter peptides and to play nonenzymatic roles in the endocrine and nervous systems65. However, previous studies have also proven that CPE plays various roles in various cancers66. Its N-terminal truncated protein has been found to be expressed in multiple cancers, regulate metastatic gene expression, and encompass different signaling pathways67. CPE has been reported to reduce aerobic glycolysis and migration in glioblastoma cells68. In addition, CPE may promote migration and invasiveness in various other cancers, such as lung cancer69, pancreatic cancer70, and cervical cancer71. In osteosarcoma, Fan et al. found that CPE-variant proteins increase invasiveness and promote epithelial-mesenchymal transition through the activation of the Wnt pathway in OS cells72. CPE also modulates immunity. Bar et al. found that it regulates immune homeostasis and inhibits inflammation73. Sanjay also reported that it enhances the innate immunity of the male reproductive tract74. Therefore, whether CPE plays a role in cancer immunity needs to be discussed further. PLB1 has seldom been studied in cancer. Its rearrangement with ALK has been found in lung cancer75, which indicates a sensitivity to targeted therapy76. Lin et al. reported that overexpressed PLB1 antigens cause high infiltration of immune cells and a favorable prognosis in glioblastoma patients77. PLB1 is commonly discussed in yeast and was found to be essential for the survival and virulence of Cryptococcus78 and Candida albicans79.

Recently, with the appearance of advanced bioinformatics algorithms, numerous gene signatures have been built from public online datasets. However, the universalism of those signatures was largely undetected or only detected in one independent cohort. Regarding the limitations of previous findings, the nine-gene signatures were summarized from the generality of several datasets. To determine the priority of our hypothesis, we included previous potential survival-related OS signatures from recent studies published in journals with a high influence and good reputation for comparison80,81,82,83,84,85,86,87,88,89,90,91,92. As we expected, after validating them in the training and validation cohorts, the nine-gene signature was the only signature showing a consistent predictive value for OS survival. Among the previous studies, the signatures by Fan et al. and Feng et al. could accurately predict the prognosis of OS in the TARGET and GSE21257 datasets. However, in the GSE39055 dataset, the lines of different groups predicted by Feng’s signature in the K‒M plot were very close, while the survival difference of groups predicted by Fan’s signature showed a reverse tendency compared with other datasets.

A GSEA of GO and KEGG pathways found that antigen processing- and presentation-related biological functions and olfactory transduction receptor signaling pathways play important roles in signature functioning. Numerous studies have proven that antigen processing- and presentation-related biological functions are key functions in the cancer immune response process and are mainly carried out by dendritic cells, macrophages and B cells93,94,95. These results indicate the important role of the three antigen-presenting cells and cancer immunity in the prognosis of OS, which is consistent with our hypothesis. Regarding the olfactory transduction receptor signaling pathway, a previous study reported that olfactory receptors provide innate and adaptive immune responses during the virus entry process96. Orecchioni et al. recently found that immune cells, such as vascular macrophages, express olfactory receptors, which induce interleukin-1β secretion, leading to inflammation97. Nevertheless, there have not been any studies that have reported any association between the olfactory transduction receptor signaling pathway and cancer immunity thus far.

Our analysis also showed the homogeneity and heterogeneity of the osteosarcoma immune microenvironment. M0 macrophages were the dominant infiltrating immune cell type in all datasets. In addition, the infiltration of M2 macrophages, which have been reported to play an immune repressive role in cancer98, was similar to that of M0 macrophages in the three datasets, except for the GSE39055 dataset. There were also several infiltrated CD8 + T cells in the TARGET and GSE16091 datasets. Other immune cell infiltration levels were quite low. The above results may partly explain the unsatisfactory response to current immunotherapy in osteosarcoma. Most of the immune microenvironment of osteosarcoma tissues consists of immune repressive M2 macrophages or nonfunctional M0 macrophages. Patients with tumor tissues infiltrated with several CD8 + T cells may receive survival benefits from immunotherapy. Given the large proportion of M0 and M2 macrophages in osteosarcoma immune cells, inducing M0 or M2 macrophages toward the M1 phenotype to promote the antitumor immune response may be a promising treatment strategy for OS patients.

There are several limitations in the present study. First, the prognostic signature was constructed from public online retrospective data with relatively small sample sizes. It may be improved by future OS studies with larger sample sizes. Second, the predictive value of the signature needs to be confirmed by future prospective studies. Third, the mechanism and function of a high infiltration of M0 macrophages in the OS microenvironment is still unclear and needs to be clarified in the future. Fourth, many prognostic factors were unavailable in the datasets used, and as such, the independence of the signature could not be fully determined. In addition, there are no wet-lab experimental data supporting the roles of signature genes and the olfactory transduction receptor signaling pathway in immune infiltration and other biological processes. Therefore, further research is needed to investigate the mechanisms.

Conclusions

In the present study, we built a new and robust nine-gene signature based on the hypothesis that the dominant infiltrated immune cells play the most important role in OS progression. The newly defined gene signature was found to be significantly associated with OS prognosis in all datasets. In addition, we proved the advantage of the signature by comparing it to previously published signatures. Antigen processing- and presentation-related biological functions and the olfactory transduction receptor signaling pathway were found to be associated with the signature risk score. The potential role and mechanism of the olfactory transduction receptor signaling pathway and M0 macrophages in OS should be evaluated in the future.