Journal of Oncology

Journal of Oncology / 2021 / Article

Research Article | Open Access

Volume 2021 |Article ID 9967954 | https://doi.org/10.1155/2021/9967954

Kang-Wen Xiao, Zhi-Bo Liu, Zi-Hang Zeng, Fei-Fei Yan, Ling-Fei Xiao, Jia-Li Li, Lin Cai, "Construction and Validation of a Macrophage-Associated Risk Model for Predicting the Prognosis of Osteosarcoma", Journal of Oncology, vol. 2021, Article ID 9967954, 18 pages, 2021. https://doi.org/10.1155/2021/9967954

Construction and Validation of a Macrophage-Associated Risk Model for Predicting the Prognosis of Osteosarcoma

Academic Editor: Raffaele Palmirotta
Received16 Mar 2021
Accepted20 May 2021
Published03 Jun 2021

Abstract

Background. Osteosarcoma is one of the most common bone tumors among children. Tumor-associated macrophages have been found to interact with tumor cells, secreting a variety of cytokines about tumor growth, metastasis, and prognosis. This study aimed to identify macrophage-associated genes (MAGs) signatures to predict the prognosis of osteosarcoma. Methods. Totally 384 MAGs were collected from GSEA software C7: immunologic signature gene sets. Differential gene expression (DGE) analysis was performed between normal bone samples and osteosarcoma samples in GSE99671. Kaplan–Meier survival analysis was performed to identify prognostic MAGs in TARGET-OS. Decision curve analysis (DCA), nomogram, receiver operating characteristic (ROC), and survival curve analysis were further used to assess our risk model. All genes from TARGET-OS were used for gene set enrichment analysis (GSEA). Immune infiltration of osteosarcoma sample was calculated using CIBERSORT and ESTIMATE packages. The independent test data set GSE21257 from gene expression omnibus (GEO) was used to validate our risk model. Results. 5 MAGs (MAP3K5, PML, WDR1, BAMBI, and GNPDA2) were screened based on protein-protein interaction (PPI), DGE, and survival analysis. A novel macrophage-associated risk model was constructed to predict a risk score based on multivariate Cox regression analysis. The high-risk group showed a worse prognosis of osteosarcoma ( < 0.001) while the low-risk group had higher immune and stromal scores. The risk score was identified as an independent prognostic factor for osteosarcoma. MAGs model for diagnosis of osteosarcoma had a better net clinical benefit based on DCA. The nomogram and ROC curve also effectively predicted the prognosis of osteosarcoma. Besides, the validation result was consistent with the result of TARGET-OS. Conclusions. A novel macrophage-associated risk score to differentiate low- and high-risk groups of osteosarcoma was constructed based on integrative bioinformatics analysis. Macrophages might affect the prognosis of osteosarcoma through macrophage differentiation pathways and bring novel sights for the progression and prognosis of osteosarcoma.

1. Introduction

Osteosarcoma, as a common malignant tumor, occurred mostly in the distal femur and proximal tibia metaphyses. Currently, the incidence of osteosarcoma worldwide was three to four per million people [1]. The main treatment for osteosarcoma consisted of chemotherapy and surgery [2]. As a highly aggressive tumor, nearly half of osteosarcoma would metastasize, and the lung was the most common metastatic site [3]. Despite a relatively low incidence of osteosarcoma, the prognosis of osteosarcoma continued to be very poor, and the cure rate after surgery was not high (17%) [4]. Hence, it was significant to explore new and effective methods to treat osteosarcoma.

The tumor microenvironment (TME) mainly consisted of tumor cells, extracellular matrix proteins, blood vessels, fibroblasts, immune cells, endothelial cells, and extracellular factors [5]. In the last few decades, TME has been paid more attention and studied in many fields, including tumor angiogenesis and metastasis [6]. Being an important part of the TME, tumor-associated macrophages (TAM) have been observed to affect the metastasis and prognosis of tumors [7]. A recent study showed that the infiltration of macrophages was related to the prognosis of breast cancer through paracrine interaction between breast cancer cells and macrophages [8]. A previous research also revealed that macrophages would promote the growth of epithelial cells with cancer-related mutations [9]. A recent study also implied that TAM prevented metastasis in high-grade osteosarcoma by collecting 132 clinical osteosarcoma samples [10]. TAM could be divided into three types based on their functional properties: M1, M2, and M0 [11]. M1 macrophages, activated by lipopolysaccharides, Th1, and other cytokines, played an important role in promoting inflammation and antimicrobial [12]. Meanwhile, M2 macrophages, mainly induced by CSF-1 and IL-10, were involved in tumor progression, wound healing, tissue repairing, and inhibition of inflammation [13]. It is currently believed that M2 macrophages promoted tumor growth and metastasis while M1 macrophages inhibited tumor formation by secreting cytokines [14]. An injection of M1 macrophages into mice with Ehrlich ascites carcinoma could improve the survival rate of the mice [15]. A former study also indicated that M1 macrophages could reduce the growth of colon cancer cells [16]. Moreover, a high proportion of M2 macrophages could lead to a poor prognosis in ovarian cancer [17]. In Zhou’s study, osteosarcoma metastasis could be prevented by all-trans retinoic acid through the prohibition of M2 polarization [18]. Although osteosarcoma and TAM have been studied in recent years, the TAM-related diagnosis and prognostic indicators of osteosarcoma were still rare. Furthermore, most of the current osteosarcoma indicators were based on the hematological study [19], and the influence of the tumor microenvironment on osteosarcoma was rarely considered. Therefore, developing a macrophage-associated risk model to predict the prognosis of osteosarcoma was urgently needed.

Therapeutically applicable research to generate effective treatments (TARGET) was a database of pediatric tumors, including acute lymphoblastic leukemia, acute myeloid leukemia, kidney tumors, neuroblastoma, and osteosarcoma. The CIBERSORT deconvolution algorithm was a machine learning method to calculate the infiltration of various immune cells based on bulk transcriptome data through linear support vector regression and highly robust to noise [20]. LM22 was a gene expression label matrix of immune cells. It consisted of 547 genes that could differentiate 22 immune cell phenotypes [21]. CIBERSORT has been widely used in predicting the proportion of immune cells in cancers. The Connectivity Map (cMap) was an instrumental online bioinformatics database that includes gene expression profiles of more than 7,000 tumor cell lines and 1,309 drugs [22]. Decision curve analysis was a statistical method to evaluate clinical prediction models, diagnostic experiments, and molecular markers [23]. Nomogram could simplify the complex prediction model in the probability estimation of the event (such as death or recurrence) based on the specific situation of each patient. Multivariate Cox regression analysis has been widely applied in clinical research [24]. Here, the clinical and transcriptome data of osteosarcoma from the TARGET database (TARGET-OS) were collected in this study. 384 MAGs were collected from GSEA software C7: immunologic signature gene sets [25, 26]. The STRING database was further used to detect MAGs with a high connection [27]. Then 5 MAGs were selected to construct the risk model through integrative bioinformatics analysis. Moreover, the prognostic nomogram was constructed to evaluate our risk model and further validated by a bootstrap test. Besides, the independent data set from the GEO database was collected to validate our model. In this study, we aimed to identify macrophage-associated gene signatures to predict the prognosis of osteosarcoma, which could help clinicians evaluate patients’ conditions and provide novel insights for osteosarcoma.

2. Materials and Methods

2.1. Sample Collection and Data Preprocessing

TARGET-OS including 88 osteosarcoma samples and clinical information was collected from the TARGET database (https://ocg.cancer.gov/programs/target). Besides, the human genome annotation GTF file was collected from the GENCODE platform (https://www.gencodegenes.org/). Moreover, the test data set GSE21257 consisting of 53 osteosarcoma samples and clinical information was collected from the GEO database (https://www.ncbi.nlm.nih.gov/geo/). TARGET-OS had 88 samples including 50 males and 37 females (1 sample with unknown gender) while GSE21257 consisted of 53 samples including 34 males and 19 females. GSE99671 had 36 samples (22 males and 14 females) including 18 normal bone samples and 18 osteosarcoma samples. The median ages (interquartile range (IQR)) of TARGET-OS and GSE21257 were 14.5 (12.2–17.4) years and 16.7 (13.6–18.7) years, respectively. GPL10295 platform was used for the GSE21257 data set while the GPL20148 platform was used for the GSE99671data set. Robust spline normalization was performed in GSE21257 while normalization of fragments per kilobase of exon model per million mapped fragments (FPKM) was performed in the TARGET-OS data set. Gene sets involving macrophage from GSEA software C7: immunologic signature were selected by searching using the keyword “macrophage.” Then, a total of 384 MAGs were collected after removing overlapped genes. These 384 MAGs were listed in Table S1. All data were normalized through the z-score method:where Zwas the standard value; xwas the specific gene expression value; µwas the mean expression value of each sample; and was the standard deviation.

The probe was a fluorescent-labeled nucleic acid complementary to a specific nucleotide sequence of the target gene. The RNA or cDNA fragment of the gene was captured by base complementary hybridization. Here, we used each probe to match genes. The maximum value of the probe was selected when the gene matched at least two probes. TARGET-OS was the training data set while GSE21257 was the test data set. Detailed information about patients with complete clinical information in TARGET-OS and GSE21257 were presented in Tables 1 and 2. Basic information of three data sets were provided in Table 3. The flowchart of this study was shown in Figure 1.


CharacteristicsRisk group value
Low-risk groupHigh-risk group

Age (IQR)15.67 (12.36–16.49)13.35 (10.08–16.23)0.173
Gender
 Female15120.311
 Male1419
Race
 White26250.384
 Black52
 Asian24
Tumor metastasis
 Metastatic480.161
 Nonmetastatic2923

IQR: interquartile range.

CharacteristicsRisk group value
Low-risk groupHigh-risk group

Age (IQR)16.67(14.59–18.59)15.08(13–18.17)0.258
Gender
 Female1090.854
 Male1717

IQR: interquartile range.

TARGET-OSGSE21257GSE99671

Osteosarcoma samples885336
Male503422
Female371914
Research objectHumanHumanHuman
Time of uploading chipAug 17,2019Mar 22,2012Nov 03, 2017
Median ages (IQR)14.5(12.2–17.4)16.7(13.6–18.7)NA

IQR: interquartile range, NA: not applicable.
2.2. Differential Gene Expression Analysis and Protein-Protein Interaction Analysis

Totally 384 MAGs were used for DGE using the limma package [28] in R software. Gene signatures of normal bone samples or osteosarcoma were identified based on DGE analysis. Benjamini–Hochberg (BH) method [29] was used here to adjust multiple hypotheses. The screening criteria for significant genes were adjusted to  < 0.05. The differentially expressed genes were further transformed into GRP files and uploaded to the cMap database (http://www.broad.mit.edu/cmap/) for small molecule drug prediction analysis. Potential therapeutic targets for osteosarcoma were selected based on enrichment score and  < 0.05. Meanwhile, 384 MAGs were used for PPI analysis in the STRING database and further visualized by Cytoscape [30]. The confidence level was 0.4, and MAGs with degrees higher than 2 were selected for further analysis.

2.3. Analysis of Infiltration of Immune Cells in Osteosarcoma Samples

The infiltration of immune cells in TARGET-OS was calculated by the CIBERSORT deconvolution algorithm. As a part of the CIBERSORT deconvolution algorithm, the machine learning method, nu-support vector regression (ν-SVR), was applied to this analysis. Different ν values including 0.25, 0.5, and 0.75 were selected to perform ν-SVR to predict the infiltration of immune cells for each sample. The solution of ν-SVR was screened based on the lowest root-mean-square error between the true and the predicted expressions. The formula of the CIBERSORT algorithm was as follows:where represented the expression level of gene i in mixed sample j, which was the sum of its expression in r immune cell type. was a label matrix (LM22; https://cibersort.stanford.edu/download.php), which represented the gene expression level of gene i in immune cells. represented the proportion of cell types in the mixed sample j. The permutations of the signature matrix were 1,000. The ESTIMATE package was also used to calculate the immune and stromal scores of each sample [31].

2.4. Survival Analysis

The survival and survminer packages were used for survival analysis in R software. The collected 384 MAGs were divided into high- and low-risk groups based on their expression in the TARGET-OS data set for survival analysis. Besides, high- and low-risk groups from the TARGET-OS and the GSE21257 data sets were used for survival analysis, respectively. High/low immune and high/low stromal score groups were also used for survival analysis. The survival rate was compared by the log-rank test.  < 0.05 indicated statistical significance.

2.5. Construction of Multivariate Cox Regression Model

Multivariate Cox regression model was constructed based on 5 differentially expressed MAGs from the TARGET-OS data set using survival package (https://cran.r-project.org/web/packages/survival/index.html) and survminer package (https://cran.r-project.org/web/packages/survminer/index.html) in R software. For each MAG, the coefficient of multivariate Cox regression was regarded as the coefficients in the risk score. Then the risk score of each sample was calculated. The formula was as follows:

where coef was the coefficient of multivariate regression analysis of MAGs. Xwas the expression of each MAG. This model was constructed to predict risk scores in osteosarcoma samples. The risk here referred to the risk of a poor prognosis in the osteosarcoma samples. The higher the risk score, the higher the probability of the poor prognosis in osteosarcoma, and vice versa. Osteosarcoma samples were divided into high- and low-risk groups by risk score. Pheatmap package (https://cran.r-project.org/web/packages/pheatmap/index.html) was used to display the expression of MAGs in the high- and low-risk groups in R software. Then risk score and other clinical characteristics were also used for multivariate Cox regression analysis to identify potential independent prognostic factors of osteosarcoma in the TARGET- OS and the GSE21257 data sets.

2.6. Decision Curve and ROC Analyses

The DCA curve was shown with the net benefit rate as the ordinate and the high-risk thresholds as the abscissa using the rmda package (https://cran.r-project.org/web/packages/rmda/index.html) in R software. The calculation formula of net benefit was as follows:

where nwas the sample size, and was the threshold probability. Here, we explored three clinical prognostic models in TARGET-OS: simple clinical data model (gender, race, age, and tumor metastasis), simple MAGs model (5 MAGs and risk score), and complex model that integrated MAGs and clinical features (gender, race, age, tumor metastasis, 5 MAGs, and risk score). The model with the greatest net benefit under different high-risk thresholds was recommended for clinical applications of osteosarcoma. Meanwhile, ROC curves to predict the 1-, 3-, and 5-year survival of TARGET-OS and GSE21257 were performed using the time ROC package in R software [32].

2.7. Construction and Internal Validation of Prognostic Nomogram

The nomogram was constructed to predict the overall survival (1 year, 3 years, and 5 years) of patients in TARGET-OS using the rms package (https://cran.r-project.org/web/packages/rms/index.html) in R software. Here, patients’ clinical characteristics such as age, gender, tumor metastasis, race, and risk score were used for the construction of a nomogram. Then internal validation of nomogram was performed, and bootstrap was set to 1,000. The discrimination of nomogram was evaluated by concordance index while calibration plots of 1-, 3-, and 5-year survival curves of TARGET-OS were performed to evaluate the prediction performance of nomogram.

2.8. Gene Set Enrichment Analysis

All genes from TARGET-OS were used for GSEA based on low- and high-risk groups. Significant macrophage-associated pathways in the osteosarcoma microenvironment were identified using GSEA software. False discovery rate (FDR) < 0.05 and  < 0.05 indicated statistical significance.

2.9. Statistical Analysis

Statistical Product and Service Solutions software (SPSS 22.0) and R 3.6.2 were used for data analysis. A chi-square test was performed for categorical data. The independent-samples Kruskal–Wallis test was performed to compare MAGs expression between normal bone samples and osteosarcoma samples. Meanwhile, the Kruskal-Walls test was also performed to compare stromal and immune scores between high- and low-risk groups, respectively. All significance levels were  < 0.05.

3. Results

3.1. Identification of 5 MAGs Related to the Prognosis of Osteosarcoma

A total of 384 MAGs were screened and used for subsequent analysis. PPI analysis of these 384 MAGs was shown in Figure 2(a). Most MAGs were well-connected to each other. Compared with macrophages in normal bone samples, TAM had different functions and played an important role in tumor progression. Therefore, DGE analysis was used to identify marker MAGs in normal bone samples and osteosarcoma. The result of DGE analysis between normal bone sample and osteosarcoma in the GSE99671 data set was shown in Figure 2(b). Gene BAMBI ( = 0.024) and gene ALOX5AP ( = 0.001) were top upregulated genes in osteosarcoma while gene WDR1 ( = 0.042), gene PML ( = 0.012), gene MAP3K5 ( = 2.60E − 05), gene GNPDA2 ( = 0.027), gene CCL5 ( = 7.37E − 08), and gene MAOA ( = 9.08E − 07) were top upregulated genes in normal bone sample. cMap analysis was also performed based on these differentially expressed MAGs, and the results were shown in Table 4. Exisulind ( = 0.031), atractyloside ( = 0.006), and doxycycline ( = 0.010) were identified as potential drugs for osteosarcoma. Immune infiltration of each sample in TARGET-OS and the correlation between immune cells were shown in Figures 2(c) and 2(d), respectively. The average proportions of M0, M1, and M2 macrophages were 43.4%, 2.35%, and 27.2%, respectively. The median proportions (IQR) of M0, M1, and M2 macrophages were 46.6% (32.2%–54.3%), 1.36% (0.420%–3.41%), and 23.4% (18.5%–36.8%), respectively. Moreover, Kaplan–Meier survival analysis of 384 MAGs in TARGET-OS was performed. Then a total of 5 intersection MAGs (WDR1, PML, MAP3K5, GNPDA2, and BAMBI) were screened based on the results of DGE, PPI, and survival analysis. The influence of different gene expression levels on the prognosis of osteosarcoma was further explored. The survival curve of each MAG in TARGET-OS was shown in Figures 3(a)–3(e). Gene WDR1 ( = 0.003), PML ( = 0.002), MAP3K5 ( < 0.001), and GNPDA2 ( = 0.030) showed a better prognosis of osteosarcoma in the high expression group. Moreover, gene BAMBI in the high expression group showed a worse prognosis ( = 0.013). Besides, the relative expression of each MAG in normal bone samples and osteosarcoma sampleswas shown in Figures 3(f)–3(j). The expression of gene WDR1 ( = 0.006), PML ( = 0.005), MAP3K5 ( = 0.010), and GNPDA2 ( = 0.009) were significantly higher in normal bone sample while in osteosarcoma samples, the expression of BAMBI ( = 0.006) was significantly higher. Then the effect of different expression levels of these 5 MAGs on gender and age was also explored. However, the results turned out that these 5 MAGs were not significantly related to gender or age. These 5 MAGs were used for subsequent analysis.


cMap nameMeannEnrichmentSpecificityPercent nonnull

Exisulind−0.5162−0.8740.031670.0287100
Chenodeoxycholic acid−0.5494−0.6980.017510.175
Atractyloside−0.4455−0.6950.005810.022760
Clorsulon−0.4094−0.6870.020770.056750
Doxycycline−0.3915−0.6640.010070.022660
Paromomycin−0.5154−0.6540.03296075
Naltrexone−0.3525−0.6430.014480.056260
Chlormezanone−0.4264−0.6410.039750.027675
Indometacin0.16580.5080.019460.015650
Flufenamic acid0.30560.5270.044790.051150
Prilocaine0.3160.5310.041770.008150
Orphenadrine0.45560.5610.026740.046166
Hydralazine0.25960.6070.0114050

3.2. Construction of a Macrophage-Associated Risk Model

The macrophage-associated risk model was constructed by 5 MAGs through multivariate Cox regression analysis in TARGET-OS. The formula of our risk model was as follows:

Then, osteosarcoma samples were divided into high- and low-risk groups by risk score, and the results were shown in Figures 4(a) and 4(b). As the risk score of patients increased, the number of deaths rose, and the survival time of patients decreased. Besides, our score could effectively predict the prognosis of osteosarcoma (concordance index = 0.797). Therefore, in order to further explore the prognostic difference between the high- and low-risk groups, clinical information of high- and low-risk groups were used as training data for survival curve analysis. The result was shown in Figure 4(c). Compared with the high-risk group, the low-risk group had a significantly improved prognosis ( < 0.001). The 5-year survival rates of high- and low-risk groups were 43.0% and 90.1%, respectively. A heat map of the expression of 5 MAGs in TARGET-OS was displayed in Figure 4(d). Besides, different stromal and immune scores among high- and low-risk groups were explored. The results turned out that the stromal and immune scores of the low-risk group were both significantly higher than the high-risk group (Figures 5(a) and 5(b)). The survival plot of each group with different immune and stromal scores was subsequently shown in Figures 5(c) and 5(d). The low-risk group continued to have a significantly better prognosis ( < 0.0001), regardless of its immune and stromal scores. Since there was a significant prognostic difference between high- and low-risk groups, our risk score was considered to predict the prognosis of osteosarcoma effectively.

3.3. Macrophage-Associated Risk Score: An Independent Prognostic Factor of Osteosarcoma

Clinical characteristics might be correlated with the prognosis of osteosarcoma. To explore whether our risk model could independently predict the prognosis, the risk score and other clinical information of TARGET-OS including gender, tumor metastatic, race, and age were used for Cox regression analysis. The results were shown in Figure 6(a). In univariate Cox regression analysis, risk score ( < 0.001) and tumor metastatic ( = 0.003) were closely related to the prognosis of osteosarcoma. Moreover, in multivariate Cox regression analysis, risk score ( < 0.001) and tumor metastatic ( = 0.001) were still related to the prognosis of osteosarcoma, indicating that tumor metastatic and our risk score could be considered independent prognostic factors of osteosarcoma.

3.4. MAG Model for Diagnosis of Osteosarcoma Had a Better Net Clinical Benefit than the Simple Clinical Model

The decision curve of three clinical prognostic models (simple clinical data model (gender, race, age, and tumor metastasis), simple MAGs model (5 MAGs and risk score), and complex model) in TARGET-OS were shown in Figures 6(b)6(d). Compared with clinical models, our MAGs model had a better net benefit for patients’ 3- and 5-year survival rate (Figures 6(b) and 6(c)). As shown in Figure 6(d), compared with clinical models, our MAGs model had the greatest net benefit to diagnose osteosarcoma metastasis. Among them, the simple clinical data model exhibited the lowest net benefit while the complex model had the highest net benefit. Therefore, a comprehensive analysis of clinical and genetic information could increase the net benefit of the model. Besides, in Figure 6(e), the area under the curve for prediction of 1-, 3-, and 5-year survival of osteosarcoma was 0.78, 0.84, and 0.84, respectively. This result also indicated that our MAGs model could accurately predict the prognosis of osteosarcoma.

3.5. Nomogram Effectively Predicted the Prognosis of Osteosarcoma

Since considering both clinical features and our MAGs model had the best net clinical benefit, the nomogram was constructed by integrating gender, age, tumor metastasis, race, and risk score. The result was displayed in Figure 7(a). As the total points became higher, the 1-, 3-, and 5-year survival rate of patients decreased. The concordance index of the nomogram was 0.842, and the calibration curve of 1-, 3-, and 5-year survival of osteosarcoma (Figures 7(b)7(d)) also illustrated our nomogram that could effectively predict the prognosis of osteosarcoma.

3.6. Validation of Our Risk Model by Independent Data Set

The independent data set GSE21257 was used for the validation of our risk model. The GSE21257 data set was divided into high- and low-risk groups based on the macrophage-associated risk model. As shown in Figures 8(a) and 8(b), compared with the low-risk group, the high-risk group had a lower survival time. The survival curve of different risk groups was shown in Figure 8(c). The 5-year survival rates of the high- and low-risk groups were 48.1% and 76.9%, respectively. The low-risk group had better clinical outcome ( = 0.040). The areas under the curve for prediction of 1-, 3-, and 5-year survival of osteosarcoma were 0.76, 0.72, and 0.73, respectively (Figure 8(d)). Moreover, the risk score and clinical information of GSE21257 including gender and age were also used for Cox regression analysis. The risk score was significantly correlated with the prognosis of osteosarcoma in univariate ( = 0.024) and multivariate ( = 0.020) regression analyses (Figure 9(a)). Therefore, our risk score could be considered an independent prognostic factor of osteosarcoma.

3.7. Low-Risk Group Was Related to Macrophage Differentiation Pathway

GSEA was also performed for all genes from TARGET-OS. The results were shown in Table S2, and important pathways are shown in Figures 9(b)9(d). Genes were enriched in GO: macrophage activation (enrichment score (ES) = 0.56,  < 0.001, and FDR < 0.001), GO: macrophage migration (ES = 0.54,  = 0.001, and FDR = 0.003), GO: macrophage differentiation (ES = 0.56,  = 0.001, and FDR = 0.004), and GO: regulation of macrophage chemotaxis (ES = 0.59,  = 0.004, and FDR = 0.009). These biological processes were considered important pathways of TAM in the osteosarcoma microenvironment. The above analysis revealed that MAGs played a significant role in osteosarcoma. Moreover, the results of characteristics of clinical information and prognosis of osteosarcoma were consistent with the result of TARGET-OS. Therefore, our risk model could differentiate different risk groups successfully, and the low-risk group was correlated with a better prognosis of osteosarcoma.

4. Discussion

As osteosarcoma is one of the most common childhood tumors in the world, its treatment and prognosis have received widespread attention. Despite the surgical treatment along with pre- and postoperative chemotherapy, the event-free 5-year survival rate remained low [33]. A previous study revealed that gene PPARG, gene IGHG3, and gene PDK1 were correlated with osteosarcoma [34]. However, this conclusion was not validated by the independent data set. In this study, DGE analysis was performed based on normal bone samples and osteosarcoma samples to identify differentially expressed MAGs. We further identified exisulind, atractyloside, and doxycycline as top potential drugs for the treatment of osteosarcoma. A recent study illustrated that exisulind could induce apoptosis in lung cancer by downregulating cyclic GMP phosphodiesterase and further improve the prognosis of lung cancer [35]. Atractyloside was also found to suppress the progression and metastasis of colon cancer in Lu’s study [36]. A previous report also showed that doxycycline could reduce the proliferation of melanoma cells by inhibiting apoptosis [36]. Considering the important role these drugs played in other tumors, they might become promising drugs for osteosarcoma treatment. The selected 5 MAGs (MAP3K5, PML, WDR1, BAMBI, and GNPDA2) were identified as prognostic-related genes for osteosarcoma, and a novel macrophage-associated risk model was constructed based on these 5 MAGs. MAP3K5, also called ASK1, was an important kinase in the process of cell apoptosis, participating in the regulation of multiple cell-signaling pathways in inflammation and tumors [37]. A recent research indicated that MAP3K5 promoted macrophage activation and migration in mice models [38]. Furthermore, knockout MAP3K5 mice were more likely to develop colon cancer [39]. These research works were also consistent with our study. MAP3K5 might be related to a better prognosis of a tumor by activating macrophages, which also needed further study to validate. The biological function of the gene PML was associated with its nuclear location, and PML was associated with cell cycle regulation and tumor suppression [40]. For instance, reduced expression of PML was found to promote multiple tumor growth, including prostate adenocarcinoma, breast carcinoma, and thyroid carcinoma [41]. A previous research reported that PML was critical to the formation of macrophages [42]. Here, we reported that high expression of PML might inhibit osteosarcoma growth, which was also consistent with these findings. WDR1 took part in inducing the disassembly of actin filaments, and dysfunction of WDR1 might cause autoinflammatory disease, which in turn activated macrophage [43]. BAMBI, a pseudoreceptor of the TGF signaling pathway, played a key role in regulating macrophage polarization [44]. A former study indicated that highly expressed BAMBI could contribute to colon cancer metastasis through Wnt/beta-catenin in mice models [45]. Besides, the expression of BAMBI was significantly higher in ovarian cancer through TGF-beta signaling [46]. BAMBI was also highly expressed in pancreatic cancer by the TGF-beta pathway [47]. Similarly, our study reported that high expression of BAMBI would result in a poor prognosis of osteosarcoma, which implied that BAMBI could be a new target for the treatment of osteosarcoma. GNPDA2 was closely related to obesity and body mass index. A recent epidemiological survey also showed that people with a high body mass index were at higher risk of cancer [48]. Besides, macrophage accumulating in adipose tissue was related to obesity [49]. Therefore, these genes were considered to be potential targets for the treatment of osteosarcoma, and further experiments are needed to verify the role of these genes in osteosarcoma.

Our study showed that the low-risk group had significantly higher immune and stromal scores, which was correlated with a better prognosis. Similarly, a previous study indicated that osteosarcoma samples with high immune scores had a better prognosis [34]. ROC curve also exhibited excellent accuracy of our risk model in TARGET-OS and GSE21257. The decision curve considered the clinical utility of specific models and focused on the net benefit of different clinical interventions. Here, the decision curve of three models (simple gene model, simple clinical model, and complex model) demonstrated that comprehensive consideration of genetic information and clinical information could obtain the greatest benefits. Besides, the nomogram also accurately predicted the prognosis of osteosarcoma. Therefore, our risk model could effectively predict the prognosis of osteosarcoma, which could be a potential clinical indicator of osteosarcoma. The GSEA results of TARGET-OS showed that GO: macrophage activation, GO: macrophage migration, GO: macrophage differentiation, and GO: regulation of macrophage chemotaxis were considered key pathways in osteosarcoma. Macrophage activation and migration were correlated with tumor growth. A recent study reported that activated TAM could promote the angiogenesis of breast cancer [50]. TAM also promoted cancer development by upregulation of LAMP2a [51]. The migration of macrophages was inhibited under hypoxia, which was conducive to tumor growth [52]. Besides, macrophage migration could accelerate tumor invasion without relying on matrix metalloproteinase [53]. Macrophage differentiation also played a significant role in the tumor microenvironment. A previous study showed that M2 macrophage was closely associated with malignant oral squamous cell carcinomas [54]. Macrophage could be shifted to M1 macrophage by phenelzine in triple-negative breast cancer [55]. Moreover, the chemokine system was found to affect macrophage polarization. A former research reported that interferon gamma could induce the activation of M1 macrophage [56]. These biological processes were closely related to TAM in tumor, which also provided a novel research perspective for the role of TAM in osteosarcoma.

Our study, however, had some limitations: (1) due to the limitation of sample size, we needed more research to support our conclusion, and (2) we also lacked prospective clinical trials to verify the performance of our model further. Therefore, we looked forward to more research on MAGs to explore novel ideas for the clinical treatment of osteosarcoma.

5. Conclusion

In general, 5 MAGs (MAP3K5, PML, WDR1, BAMBI, and GNPDA2) correlated with the prognosis of osteosarcoma were screened for the construction of the risk model. A novel macrophage-associated risk score to differentiate low- and high-risk groups of osteosarcoma was constructed based on multiple bioinformatics analyses. The high score indicated the poor prognosis of osteosarcoma while the low score indicated the better prognosis of osteosarcoma. Besides, our risk score was validated by the independent data set successfully, and nomogram effectively predicted the prognosis of osteosarcoma.

Abbreviations

MAG:Macrophage-associated gene
GSEA:Gene set enrichment analysis
DGE:Differential gene expression
PPI:Protein-protein interaction
cMap:Connectivity Map
ROC:Receiver operating characteristic
GEO:Gene expression omnibus
TME:Tumor microenvironment
TAM:Tumor-associated macrophage
TARGET:Therapeutically applicable research to generate effective treatment
FDR:False discovery rate
DCA:Decision curve analysis
ES:Enrichment score.

Data Availability

The data sets (TARGET-OS, GSE21257, and GSE99671) analyzed during the current study are available in the GEO database (https://www.ncbi.nlm.nih.gov/) and TARGET database (https://ocg.cancer.gov/programs/target).

Conflicts of Interest

The authors declare that there are no conflicts of interest.

Authors’ Contributions

L.C., K.-W.X., and Z.-B.L. conceived and designed the study. Z.-B.L. and F.-F.Y. searched databases and analyzed data. K.-W.X., J.-L.L., Z.-H.Z., and L.-F.X. prepared tables and figures. L.C., K.-W.X., and Z.-H. Z. wrote the manuscript. L.C. supervised the study. Kang-Wen Xiao and Zhi-Bo Liu contributed equally to this work.

Acknowledgments

This work was supported by the Health Care of Yellow Crane Talent Plan (Project no. 17).

Supplementary Materials

Table S1. Full list of 384 macrophage-associated genes. Table S2. Enrichment analysis in TARGET-OS. (Supplementary Materials)

References

  1. M. Kansara, M. W. Teng, M. J. Smyth, and D. M. Thomas, “Translational biology of osteosarcoma,” Nature Reviews Cancer, vol. 14, no. 11, pp. 722–735, 2014. View at: Publisher Site | Google Scholar
  2. L. Mirabello, R. J. Troisi, and S. A. Savage, “Osteosarcoma incidence and survival rates from 1973 to 2004,” Cancer, vol. 115, no. 7, pp. 1531–1543, 2009. View at: Publisher Site | Google Scholar
  3. A. J. Chou, E. S. Kleinerman, M. D. Krailo et al., “Addition of muramyl tripeptide to chemotherapy for patients with newly diagnosed metastatic osteosarcoma,” Cancer, vol. 115, no. 22, pp. 5339–5348, 2009. View at: Publisher Site | Google Scholar
  4. M. P. Link, A. M. Goorin, A. W. Miser et al., “The effect OF adjuvant chemotherapy ON relapse-free survival IN patients with osteosarcoma OF the extremity,” New England Journal of Medicine, vol. 314, no. 25, pp. 1600–1606, 1986. View at: Publisher Site | Google Scholar
  5. D. Hanahan and R. A. Weinberg, “Hallmarks of cancer: the next generation,” Cell, vol. 144, no. 5, pp. 646–674, 2011. View at: Publisher Site | Google Scholar
  6. D. Hanahan and L. M. Coussens, “Accessories to the crime: functions of cells recruited to the tumor microenvironment,” Cancer Cell, vol. 21, no. 3, pp. 309–322, 2012. View at: Publisher Site | Google Scholar
  7. J. W. Pollard, “Tumour-educated macrophages promote tumour progression and metastasis,” Nature Reviews Cancer, vol. 4, no. 1, pp. 71–78, 2004. View at: Publisher Site | Google Scholar
  8. J. B. Wyckoff, Y. Wang, E. Y. Lin et al., “Direct visualization of macrophage-assisted tumor cell intravasation in mammary tumors,” Cancer Research, vol. 67, no. 6, pp. 2649–2656, 2007. View at: Publisher Site | Google Scholar
  9. B.-Z. Qian and J. W. Pollard, “Macrophage diversity enhances tumor progression and metastasis,” Cell, vol. 141, no. 1, pp. 39–51, 2010. View at: Publisher Site | Google Scholar
  10. E. P. Buddingh, M. L. Kuijjer, R. A. J. Duim et al., “Tumor-infiltrating macrophages are associated with metastasis suppression in high-grade osteosarcoma: a rationale for treatment with macrophage activating agents,” Clinical Cancer Research, vol. 17, no. 8, pp. 2110–2119, 2011. View at: Publisher Site | Google Scholar
  11. C. Mills, “M1 and M2 macrophages: oracles of Health and disease,” Critical Reviews in Immunology, vol. 32, no. 6, pp. 463–488, 2012. View at: Publisher Site | Google Scholar
  12. A. J. Fleetwood, T. Lawrence, J. A. Hamilton, and A. D. Cook, “Granulocyte-macrophage colony-stimulating factor (CSF) and macrophage CSF-dependent macrophage phenotypes display differences in cytokine profiles and transcription factor activities: implications for CSF blockade in inflammation,” The Journal of Immunology, vol. 178, no. 8, pp. 5245–5252, 2007. View at: Publisher Site | Google Scholar
  13. A. Sica, P. Larghi, A. Mancino et al., “Macrophage polarization in tumour progression,” Seminars in Cancer Biology, vol. 18, no. 5, pp. 349–355, 2008. View at: Publisher Site | Google Scholar
  14. A. Mantovani and A. Sica, “Macrophages, innate immunity and cancer: balance, tolerance, and diversity,” Current Opinion in Immunology, vol. 22, no. 2, pp. 231–237, 2010. View at: Publisher Site | Google Scholar
  15. J. Jackute, M. Zemaitis, D. Pranys et al., “Distribution of M1 and M2 macrophages in tumor islets and stroma in relation to prognosis of non-small cell lung cancer,” BMC Immunology, vol. 19, 2018. View at: Publisher Site | Google Scholar
  16. S. V. Kalish, S. V. Lyamina, E. A. Usanova, E. B. Manukhina, N. P. Larionov, and I. Y. Malyshev, “Macrophages reprogrammed in vitro towards the M1 phenotype and activated with LPS extend lifespan of mice with Ehrlich ascites carcinoma,” Medical Science Monitor Basic Research, vol. 21, pp. 226–234, 2015. View at: Publisher Site | Google Scholar
  17. C. Le Page, A. Marineau, P. K. Bonza et al., “BTN3A2 expression in epithelial ovarian cancer is associated with higher tumor infiltrating T cells and a better prognosis,” PLoS One, vol. 7, no. 6, Article ID e38541, 2012. View at: Publisher Site | Google Scholar
  18. X.-j. Shao, S.-f. Xiang, Y.-q. Chen et al., “Inhibition of M2-like macrophages by all-trans retinoic acid prevents cancer initiation and stemness in osteosarcoma cells,” Acta Pharmacologica Sinica, vol. 40, no. 10, pp. 1343–1350, 2019. View at: Publisher Site | Google Scholar
  19. Q. Zhang, G. Dong, F. Wang, and W. Ding, “Correlation between the changes of serum COX 2, APE1, VEGF, TGF-beta and TSGF levels and prognosis in patients with osteosarcoma before and after treatment,” Journal of Cancer Research and Therapeutics, vol. 16, no. 2, pp. 335–342, 2020. View at: Publisher Site | Google Scholar
  20. A. M. Newman, C. L. Liu, M. R. Green et al., “Robust enumeration of cell subsets from tissue expression profiles,” Nature Methods, vol. 12, no. 5, pp. 453–457, 2015. View at: Publisher Site | Google Scholar
  21. B. Chen, M. S. Khodadoust, C. L. Liu, A. M. Newman, and A. A. Alizadeh, “Profiling tumor infiltrating immune cells with CIBERSORT,” Methods in Molecular Biology, vol. 1711, pp. 243–259, 2018. View at: Publisher Site | Google Scholar
  22. J. Lamb, E. D. Crawford, D. Peck et al., “The connectivity map: using gene-expression signatures to connect small molecules, genes, and disease,” Science, vol. 313, no. 5795, pp. 1929–1935, 2006. View at: Publisher Site | Google Scholar
  23. A. J. Vickers and E. B. Elkin, “Decision curve analysis: a novel method for evaluating prediction models,” Medical Decision Making, vol. 26, no. 6, pp. 565–574, 2006. View at: Publisher Site | Google Scholar
  24. F. E. Harrell, K. L. Lee, and D. B. Mark, “Multivariable prognostic models: issues in developing models, evaluating assumptions and adequacy, and measuring and reducing errors,” Statistics in Medicine, vol. 15, no. 4, pp. 361–387, 1996. View at: Publisher Site | Google Scholar
  25. A. Subramanian, P. Tamayo, V. K. Mootha et al., “Gene set enrichment analysis: a knowledge-based approach for interpreting genome-wide expression profiles,” Proceedings of the National Academy of Sciences, vol. 102, no. 43, pp. 15545–15550, 2005. View at: Publisher Site | Google Scholar
  26. V. K. Mootha, C. M. Lindgren, K.-F. Eriksson et al., “PGC-1α-responsive genes involved in oxidative phosphorylation are coordinately downregulated in human diabetes,” Nature Genetics, vol. 34, no. 3, pp. 267–273, 2003. View at: Publisher Site | Google Scholar
  27. D. Szklarczyk, A. Franceschini, S. Wyder et al., “STRING v10: protein-protein interaction networks, integrated over the tree of life,” Nucleic Acids Research, vol. 43, no. D1, pp. D447–D452, 2015. View at: Publisher Site | Google Scholar
  28. M. E. Ritchie, B. Phipson, D. Wu et al., “Limma powers differential expression analyses for RNA-sequencing and microarray studies,” Nucleic Acids Res, vol. 43, 2015. View at: Publisher Site | Google Scholar
  29. J. D. Storey, “A direct approach to false discovery rates,” Journal of the Royal Statistical Society: Series B (Statistical Methodology), vol. 64, no. 3, pp. 479–498, 2002. View at: Publisher Site | Google Scholar
  30. P. Shannon, A. Markiel, O. Ozier et al., “Cytoscape: a software environment for integrated models of biomolecular interaction networks,” Genome Research, vol. 13, no. 11, pp. 2498–2504, 2003. View at: Publisher Site | Google Scholar
  31. K. Yoshihara, M. Shahmoradgoli, E. Martinez et al., “Inferring tumour purity and stromal and immune cell admixture from expression data,” Nature Communications, vol. 4, 2013. View at: Publisher Site | Google Scholar
  32. P. J. Heagerty, T. Lumley, and M. S. Pepe, “Time-dependent ROC curves for censored survival data and a diagnostic marker,” Biometrics, vol. 56, no. 2, pp. 337–344, 2000. View at: Publisher Site | Google Scholar
  33. L. Kager, G. Tamamyan, and S. Bielack, “Novel insights and therapeutic interventions for pediatric osteosarcoma,” Future Oncology, vol. 13, no. 4, pp. 357–368, 2017. View at: Publisher Site | Google Scholar
  34. C. Zhang, J.-H. Zheng, Z.-H. Lin et al., “Profiles of immune cell infiltration and immune-related genes in the tumor microenvironment of osteosarcoma,” Aging, vol. 12, no. 4, pp. 3486–3501, 2020. View at: Publisher Site | Google Scholar
  35. C. M. Whitehead, K. A. Earle, J. Fetter et al., “Exisulind-induced apoptosis in a non-small cell lung cancer orthotopic lung tumor model augments docetaxel treatment and contributes to increased survival,” Molecular Cancer Therapeutics, vol. 2, no. 5, pp. 479–488, 2003. View at: Google Scholar
  36. L. Qi, F. Song, Y. Han, Y. Zhang, and Y. Ding, “Atractyloside targets cancer-associated fibroblasts and inhibits the metastasis of colon cancer,” Annals of Translational Medicine, vol. 8, no. 21, p. 1443, 2020. View at: Publisher Site | Google Scholar
  37. M. Yin, H. J. Zhou, J. Zhang et al., “ASK1-dependent endothelial cell activation is critical in ovarian cancer growth and metastasis,” JCI Insight, vol. 2, 2017. View at: Publisher Site | Google Scholar
  38. N. Osaka, T. Takahashi, S. Murakami et al., “ASK1-dependent recruitment and activation of macrophages induce hair growth in skin wounds,” Journal of Cell Biology, vol. 176, no. 7, pp. 903–909, 2007. View at: Publisher Site | Google Scholar
  39. Y. Hayakawa, Y. Hirata, H. Nakagawa et al., “Apoptosis signal-regulating kinase 1 regulates colitis and colitis-associated tumorigenesis by the innate immune responses,” Gastroenterology, vol. 138, no. 3, pp. 1055–1067, 2010. View at: Publisher Site | Google Scholar
  40. S. Nisole, M. A. Maroui, X. H. Mascle, M. Aubry, and M. K. Chelbi-Alix, “Differential roles of PML isoforms,” Frontiers in Oncology, vol. 3, p. 125, 2013. View at: Publisher Site | Google Scholar
  41. C. Gurrieri, P. Capodieci, R. Bernardi et al., “Loss of the tumor suppressor PML in human cancers of multiple histologic origins,” JNCI Journal of the National Cancer Institute, vol. 96, no. 4, pp. 269–279, 2004. View at: Publisher Site | Google Scholar
  42. Y. Khalfin-Rabinovich, A. Weinstein, and B.-Z. Levi, “PML is a key component for the differentiation of myeloid progenitor cells to macrophages,” International Immunology, vol. 23, no. 4, pp. 287–296, 2011. View at: Publisher Site | Google Scholar
  43. M. L. Kim, J. J. Chae, Y. H. Park et al., “Aberrant actin depolymerization triggers the pyrin inflammasome and autoinflammatory disease that is dependent on IL-18, not IL-1β,” Journal of Experimental Medicine, vol. 212, no. 6, pp. 927–938, 2015. View at: Publisher Site | Google Scholar
  44. D. Wang, X. Chen, and R. Zhang, “BAMBI promotes macrophage proliferation and differentiation in gliomas,” Molecular Medicine Reports, vol. 17, pp. 3960–3966, 2018. View at: Publisher Site | Google Scholar
  45. N. Togo, S. Ohwada, S. Sakurai et al., “Prognostic significance of BMP and activin membranebound inhibitor in colorectal cancer,” World Journal of Gastroenterology, vol. 14, no. 31, pp. 4880–4888, 2008. View at: Publisher Site | Google Scholar
  46. D. Pils, M. Wittinger, M. Petz et al., “BAMBI is overexpressed in ovarian cancer and co-translocates with Smads into the nucleus upon TGF-ß treatment,” Gynecologic Oncology, vol. 117, no. 2, pp. 189–197, 2010. View at: Publisher Site | Google Scholar
  47. J. A. Brosnan, S. Yachida, and C. A. Iacobuzio-Donahue, “BAMBI Is overexpressed in metastatic pancreatic cancers with genetically Intact TGF-beta pathways: a potential mechanism to escape TGF-beta signaling during metastasis formation,” Cancer Research, vol. 718, 2011. View at: Google Scholar
  48. M. Benn, A. Tybjærg-Hansen, G. D. Smith, and B. G. Nordestgaard, “High body mass index and cancer risk-a Mendelian randomisation study,” European Journal of Epidemiology, vol. 31, no. 9, pp. 879–892, 2016. View at: Publisher Site | Google Scholar
  49. S. Subramanian and A. Chait, “The effect of dietary cholesterol on macrophage accumulation in adipose tissue: implications for systemic inflammation and atherosclerosis,” Current Opinion in Lipidology, vol. 20, no. 1, pp. 39–44, 2009. View at: Publisher Site | Google Scholar
  50. J. S. Lewis, R. J. Landers, J. C. E. Underwood, A. L. Harris, and C. E. Lewis, “Expression of vascular endothelial growth factor by macrophages is up-regulated in poorly vascularized areas of breast carcinomas,” The Journal of Pathology, vol. 192, no. 2, pp. 150–158, 2000. View at: Publisher Site | Google Scholar
  51. R. Wang, Y. Liu, L. Liu et al., “Tumor cells induce LAMP2a expression in tumor-associated macrophage for cancer progression,” EBioMedicine, vol. 40, pp. 118–134, 2019. View at: Publisher Site | Google Scholar
  52. M. J. Grimshaw and F. R. Balkwill, “Inhibition of monocyte and macrophage chemotaxis by hypoxia and inflammation - a potential mechanism,” European Journal of Immunology, vol. 31, no. 2, pp. 480–489, 2001. View at: Publisher Site | Google Scholar
  53. R. Guiet, E. Van Goethem, C. Cougoule et al., “The process of macrophage migration promotes matrix metalloproteinase-independent invasion by tumor cells,” The Journal of Immunology, vol. 187, no. 7, pp. 3806–3814, 2011. View at: Publisher Site | Google Scholar
  54. F. Wehrhan, M. Buettner-Herold, P. Hyckel et al., “Increased malignancy of oral squamous cell carcinomas (oscc) is associated with macrophage polarization in regional lymph nodes - an immunohistochemical study,” BMC Cancer, vol. 14, 2014. View at: Publisher Site | Google Scholar
  55. A. H. Y. Tan, W. Tu, R. McCuaig et al., “Lysine-specific histone demethylase 1A regulates macrophage polarization and checkpoint molecules in the tumor microenvironment of triple-negative breast cancer,” Frontiers in Immunology, vol. 10, 2019. View at: Publisher Site | Google Scholar
  56. L. Zhu, Q. Zhao, T. Yang, W. Ding, and Y. Zhao, “Cellular metabolism and macrophage functional polarization,” International Reviews of Immunology, vol. 34, no. 1, pp. 82–100, 2015. View at: Publisher Site | Google Scholar

Copyright © 2021 Kang-Wen Xiao et al. This is an open access article distributed under the Creative Commons Attribution License, which permits unrestricted use, distribution, and reproduction in any medium, provided the original work is properly cited.

Related articles

No related content is available yet for this article.
 PDF Download Citation Citation
 Download other formatsMore
 Order printed copiesOrder
Views287
Downloads391
Citations

Related articles

No related content is available yet for this article.

Article of the Year Award: Outstanding research contributions of 2021, as selected by our Chief Editors. Read the winning articles.