Journal of Oncology
Volume 2021 (2021), Article ID 9967954, 18 pages
https://doi.org/10.1155/2021/9967954

Research Article

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

Kang-Wen Xiao,¹ Zhi-Bo Liu,¹ Zi-Hang Zeng,² Fei-Fei Yan,¹ Ling-Fei Xiao,¹ Jia-Li Li,² and Lin Cai¹

¹Department of Orthopedics, Zhongnan Hospital of Wuhan University, Wuhan, Hubei430071, China
²Department of Radiation and Medical Oncology, Zhongnan Hospital of Wuhan University, Wuhan, Hubei430071, China

Correspondence should be addressed to Lin Cai; orthopedics@whu.edu.cn

Received 16 March 2021; Accepted 20 May 2021; Published 3 June 2021

Academic Editor: Raffaele Palmirotta

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.

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.

Table 1: Detailed clinical information about the TARGET-OS data set.

Characteristics
Risk group
value

Low-risk group
High-risk group

Age (IQR)
15.67 (12.36–16.49)
13.35 (10.08–16.23)
0.173

Gender

    Female
15
12
0.311

    Male
14
19

Race

    White
26
25
0.384

    Black
5
2

    Asian
2
4

Tumor metastasis

    Metastatic
4
8
0.161

    Nonmetastatic
29
23

IQR: interquartile range.

Table 2: Detailed clinical information about the GSE21257 data set.

Characteristics
Risk group
value

Low-risk group
High-risk group

Age (IQR)
16.67(14.59–18.59)
15.08(13–18.17)
0.258

Gender

    Female
10
9
0.854

    Male
17
17

IQR: interquartile range.

Table 3: Detailed information about the GSE21257, the GSE99671, and the TARGET-OS data sets.

TARGET-OS
GSE21257
GSE99671

Osteosarcoma samples
88
53
36

Male
50
34
22

Female
37
19
14

Research object
Human
Human
Human

Time of uploading chip
Aug 17,2019
Mar 22,2012
Nov 03, 2017

Median ages (IQR)
14.5(12.2–17.4)
16.7(13.6–18.7)
NA

IQR: interquartile range, NA: not applicable.

Figure 1: Flowchart of this study.

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.

(a)

(b)

(c)

(d)

Figure 2: Identification of prognostic MAGs through PPI and DGE and immune infiltration of TARGET-OS: (a) PPI analysis of 384 MAGs, (b) differential gene expression analysis between osteosarcoma samples and normal bone sample, (c) correlation plot of each immune cell in TARGET-OS, and (d) immune infiltration of each sample in TARGET-OS.

Table 4: Potential drugs for the treatment of osteosarcoma by cMap analysis.

cMap name
Mean
n
Enrichment

Specificity
Percent nonnull

Exisulind
−0.516
2
−0.874
0.03167
0.0287
100

Chenodeoxycholic acid
−0.549
4
−0.698
0.01751
0.1
75

Atractyloside
−0.445
5
−0.695
0.00581
0.0227
60

Clorsulon
−0.409
4
−0.687
0.02077
0.0567
50

Doxycycline
−0.391
5
−0.664
0.01007
0.0226
60

Paromomycin
−0.515
4
−0.654
0.03296
0
75

Naltrexone
−0.352
5
−0.643
0.01448
0.0562
60

Chlormezanone
−0.426
4
−0.641
0.03975
0.0276
75

Indometacin
0.165
8
0.508
0.01946
0.0156
50

Flufenamic acid
0.305
6
0.527
0.04479
0.0511
50

Prilocaine
0.31
6
0.531
0.04177
0.0081
50

Orphenadrine
0.455
6
0.561
0.02674
0.0461
66

Hydralazine
0.259
6
0.607
0.0114
0
50

(a)

(b)

(c)

(d)

(e)

(f)

(g)

(h)

(i)

(j)

Figure 3: The results of survival curve for 5 MAGs (BAMBI, PML, GNPDA2, WDR1, and MAP3K5) and the expression of each MAG in a normal bone sample and osteosarcoma: (a)–(e) survival curve of 5 MAGs (BAMBI, PML, GNPDA2, WDR1, and MAP3K5) by Kaplan–Meier method and (f)–(j) the boxplots of expression of each MAG in normal bone sample and osteosarcoma.  < 0.05 and  < 0.01.

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.

(a)

(b)

(c)

(d)

Figure 4: Construction of macrophage-associated risk model and survival curve of high-/low-risk groups in TARGET-OS: (a) TARGET-OS was divided into high- and low-risk groups using the median risk score as the cutoff value, (b) the relationship between risk score and survival time and status of patients, (c) the survival curve of high- and low-risk groups in TARGET-OS, and (d) the heat map between expression of 5 MAGs and osteosarcoma samples.

(a)

(b)

(c)

(d)

Figure 5: Comparison of stromal and immune scores among high- and low-risk groups and survival analysis: (a) comparison of stromal score among high- and low-risk groups, (b) comparison of immune score among high- and low-risk groups, (c) survival analysis of different stromal scores among high- and low-risk groups (L: low, H: high, LSS: low stromal score, and HSS: high stromal score), and (d) survival analysis of different immune scores among high- and low-risk groups (L: low, H: high, LIS: low immune score, and HIS: high immune score). < 0.05and  < 0.01.

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.

(a)

(b)

(c)

(d)

(e)

Figure 6: Identification of risk score as an independent prognostic factor of osteosarcoma, decision curve analysis of three models (simple gene model, simple clinical model, and complex model), and ROC curve for predicting the prognosis of osteosarcoma in 1, 3, and 5 years: (a) identification of risk score as an independent prognostic factor of osteosarcoma by univariate and multivariate Cox regression analyses, (b) the decision curve of the net benefit of the 2 models for the 3-year survival rate (simple gene and simple clinical models), (c) the decision curve of the net benefit of the 3 models for the 5-year survival rate (simple gene, simple clinical, and complex models); (d) the decision curve of the net benefit of the 3 models for the diagnosis of osteosarcoma metastasis (simple gene, simple clinical, and complex models), and (e) ROC curve to predict the prognosis of osteosarcoma in 1, 3, and 5 years.

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.

(a)

(b)

(c)

(d)

Figure 7: Construction and internal validation of nomogram to predict the prognosis of osteosarcoma: (a) construction of the nomogram by collecting age, gender, tumor metastasis, race, and risk score; (b) 1-year survival calibration curve; (c) 3-year survival calibration curve; and (d) 5-year survival calibration curve.

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.

(a)

(b)

(c)

(d)

Figure 8: Validation of macrophage-associated risk model and survival curve of high-/low-risk groups in GSE21257: (a) GSE21257 was divided into high- and low-risk groups using the median risk score as the cutoff value, (b) the relationship between risk score and survival time and status of patients, (c) the survival curve of high- and low-risk groups in GSE21257, and (d) ROC curve to predict the prognosis of osteosarcoma in 1, 3, and 5 years.

(a)

(b)

(c)

(d)

Figure 9: Validation of risk score as an independent prognostic factor of osteosarcoma in GSE21257 and important pathways identified by GSEA: (a) validation of risk score as an independent prognostic factor of osteosarcoma in GSE21257 by univariate and multivariate Cox regression analyses, (b) the important enrichment pathway: GO – macrophage migration, (c) the important enrichment pathway: GO-macrophage differentiation, and (d) the important enrichment pathways: GO-macrophage activation.

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

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.

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.

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.

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.

D. Hanahan and R. A. Weinberg, “Hallmarks of cancer: the next generation,” Cell, vol. 144, no. 5, pp. 646–674, 2011.

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.

J. W. Pollard, “Tumour-educated macrophages promote tumour progression and metastasis,” Nature Reviews Cancer, vol. 4, no. 1, pp. 71–78, 2004.

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.

B.-Z. Qian and J. W. Pollard, “Macrophage diversity enhances tumor progression and metastasis,” Cell, vol. 141, no. 1, pp. 39–51, 2010.

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.

C. Mills, “M1 and M2 macrophages: oracles of Health and disease,” Critical Reviews in Immunology, vol. 32, no. 6, pp. 463–488, 2012.

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.

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.

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.

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.

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.

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.

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.

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.

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.

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.

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.

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.

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.

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.

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.

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.

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.

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.

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.

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.

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.

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.

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.

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.

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.

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.

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.

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.

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.

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.

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.

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.

D. Wang, X. Chen, and R. Zhang, “BAMBI promotes macrophage proliferation and differentiation in gliomas,” Molecular Medicine Reports, vol. 17, pp. 3960–3966, 2018.

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.

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.

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.

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.

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.

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.

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.

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.

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.

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.

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.

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.