Revealing the link between macrophage in microenvironment of osteosarcoma and poor prognosis by utilizing the Integrated analysis.

OBJECTIVES: Osteosarcoma (OS) is the most common type of primary malignant bone tumor, The effect of tumor microenvironment components on OS oncogenesis remains unknown.
METHODS: To investigate the function of immune cells in osteosarcoma, we provided a text-based GMT (Gene Matrix Transposed) file in which each line defines one of lm22 with their markers. We used STRING to draw DEG's PPI network and selected hub genes and modules. Then, survival analysis was conducted to hub genes. We identified 10,390 common genes, and identified 218 DEGs based on the combined t-value and Z scores.
RESULTS: The KEGG and GSEA enrichment analysis showed that macrophages are significantly activated in osteosarcoma. PPI network analysis revealed that hub gene CD163 molecule. We found that the expression of CD163 was negatively associated with the OS of osteosarcoma patients. These results suggest that macrophages are a risk factor in patients with osteosarcoma.
CONCLUSIONS: This study has systematically validated results of the studies carried out previously and filled up the gap in the field of OS on large-scaled meta-analysis. In addition, for the hub gene (CD163) and the macrophage cell capable of being used as a novel biomarker in promoting early diagnosis and development of therapeutic approaches.

Introduction

As the most common type of primary malignant bone tumor, osteosarcoma (OS) has a poor prognosis after treatment with a combination of surgery and chemotherapy[1]. Local pain, together with local swelling and limitation in joint motion, constitutes the typical signs and symptoms of osteosarcoma and causes great pain to the patients in their daily lives. If untreated, osteosarcoma would progress relentlessly with both local and systemic pain, and result in death within a few months. Even though modern multimodality treatment could significantly improve both tumor resectability and long-term prognosis of the patients, there are still 25-35% of the patients with initially non-metastatic disease tending to develop metastasis subsequently. Therefore, to avoid chemo-resistance, reduce lung metastasis, improving clinical outcomes, and understanding the mechanisms underlying disease progression, the interaction between the different cellular players of the tumor microenvironment (TME), along with the extracellular matrix (ECM), recapitulation them is one strategy[2].

At present, DNA microarray technology has attracted great interest among researchers. Known as “gene chip” technology also, it integrates molecular genetics with computer science on a great scale. The technology can demonstrate the simultaneous expression in details of the entire genome quickly. Since high-throughput technology has been used in many fields, detecting abnormal genomic changes in OS patients via OS chips by detecting the expression levels throughout the genome is a useful way. Recently, many researchers have used this technology to increase knowledge on changes in OS in terms of cells and molecules in a more comprehensive way. In 2014, Yang et al.[3] carried out a meta-analysis of the OS microarray data to understand the underlying OS mechanism in a better way. As for this study, the data acquired from different microarray platforms and the results of OS tissue samples and cell lines were selected as the data source. It is demonstrated in this study that either the “extracellular matrix (ECM)-receptor interaction” or the “cell cycle” was a highly enriched KEGG pathway, and PTBP2, RGS4, and FXYD6, and several other hub genes were also identified in this study[3]. Either the overexpression of c-fos (FOS) or the runt-related transcription factor 2 (RUNX2) might exert a specific effect in OS; and in the patients with OS, expression of RUNX2, in particular, might act as a marker of their chemotherapy failure[4,5].

Although these studies show a list of DEGs, there are often inconsistencies among the studies because of the small sample size and the limitations of different results acquired (by other laboratory procedures, microarray platforms, and analytical technologies) from different groups. Recent studies have shown that integrating the expression data of genes from multiple sources systematically, known as a meta-analysis, could improve the statistical ability to detect DEGs and assess heterogeneity and possibly lead to more robust predictions reproducible and more accurate.

Therefore, we integrated large volumes of DNA microarray data sets for meta-analysis. We collected 6 data sets in the GEO database and performed a meta-analysis using the expression matrix. 218 differential genes were obtained through the integration analysis of t-test and z-score. Through the analysis of KEGG and GO pathways, it was found that macrophage-related pathways were activated. However, the macrophages marker (CD163) is highly expressed in osteosarcoma and is associated with prognosis. Also, through analyzing EPIC of the four RNA-seq data sets screened from the database, we found that macrophages increased in tumors compared with para-carcinoma cells. Combined with the above analysis, we found increased macrophages in osteosarcoma and confirmed that the tumor’s immune microenvironment correlated.

Material and methods

Selection of microarray datasets

Referring to the guidelines for Preferred Reporting Items for Systematic Reviews and Meta-Analyses (PRISMA), which was published in 2009, we searched the in the National Center for Biotechnology Information (NCBI, website: http://www.ncbi.nlm.nih.gov/geo/) for databases of Gene Expression Omnibus (GEO) in a comprehensive way. We used the keyword “osteosarcoma” in our search. It’s required that the datasets should meet the criteria as follows: 1) the study organisms must be Homo sapiens; 2) the datasets must contain both samples of osteosarcoma patients and samples of corresponding normal tissues. The studies would be excluded in the following circumstances: 1) they were studies on cell lines; 2) they were free of traceability in literature; 3) they were studies on DNA methylation; 4) they were studies based on miRNA; and 5) they lacked both cases and controls. In fact, two independent analysts have searched the original studies and selected the data required, and any disagreement between them was finally solved by a third analyst to which they referred for consultation.

Meta-analysis of multiple microarray datasets

It was from GEO’s database that files (.CEL) on microarray datasets meeting the criteria for inclusion were downloaded. Using R statistical software (website: http://www.r-project.org/) to integrate the genes from datasets, common genes were obtained. Then, using R statistical software, a meta-analysis of expression profiles of common genes was performed by a combination of p-values and Z scores. The combination of p-values (pval_test) has both test-statistic and p-value included. Common genes were performed a meta-analysis by packages of MAMA, mataMA, affyPLM, CLL, and RankProd. In using R statistical software to perform two meta-analyses, we used the combination meta-analysis of p-values (with the threshold being absolute value greater than 5) and Z-scores (with the threshold being absolute value greater than 4) to serve as the cut-off criteria and identified a list of (upregulated or downregulated) DEGs.

Analysis of GO annotations and KEGG pathway enrichment

It is of vital importance to identify the biological characteristics of DEGs. Based on the meta-analysis results, enrichment analysis was used to evaluate the most significant DEGs. Then, WEB-based Gene Set Analysis Toolkit (website: http://www.webgestalt.org/option.php) was applied to perform analysis of Kyoto Encyclopedia of Genes and Genomes (KEGG) pathway enrichment and annotation of gene ontology (GO), to identify the DEGs that were of the most significance.

GSEA pathway analysis and EPIC

To investigate the function of immune cells in osteosarcoma, we provided a text-based GMT (Gene Matrix Transposed) file, in which one of lm22 was defined for each line with their markers[6]. Then pathway analysis of Gene Set Enrichment Analysis (GSEA) related to genes ranked by t-values in the meta-analysis was performed. Lm22 were as follows: the B_cells_naive, the B_cells_memory, the Plasma_cells, the T_cells_CD8, the T_cells_CD4_naive, the T_cells_CD4_memory_resting, the T_cells_CD4_memory_activated, the T_cells_follicular_helper, the T_cells_regulatory_(Tregs), the T_cells_gamma_delta, the NK_cells_resting, the NK_cells_activated, the Monocytes, the Macrophages_M0, the Macrophages_M1, the Macrophages_M2, the Dendritic_cells_resting, the Dendritic_cells_activated, the Mast_cells_resting, the Mast_cells_activated, the Eosinophils and the Neutrophils. To investigate the proportion of immune cell components in osteosarcoma, we put four RNA-Seq datasets from GEO to a tool named EIPC[7]. The tool is known as “Estimate the Proportion of Immune and Cancer cells”, or EPIC for short. This tool compares the level of gene expression in a tumor and a library of the gene expression profiles from specific cell types as can be found in tumors and uses the comparison result to predict the amount of each type of cells that exist. EPIC’s predictions have been confirmed to be accurate by experimental measurements of several human tumors.

Construction of PPI network

The database of Search Tool for the Retrieval of Interacting Genes (known as STRING) (website: http://string-db.org) provides data on protein-protein interactions (known as PPIs)[8]. With the confidence score ≥0.7 as the criterion for significance cut-off, we used STRING to chart a PPI network of DEGs, aiming to have a profound understanding and prediction of the identified DEGs’ cellular functions and biological behaviors. Also, Cytoscape software was used for visualization of the PPI network[9].

Selection of hub genes and modules

PPI network was conducted calculation in terms of its degree, closeness, and betweenness using Cytoscape software. The degree of a node refers to the average number of edges (or interactions) occurring on it[10]. In accordance with the degree of a node, identification of the hub genes was made. For selecting the clustering modules of PPI network with the most significance in Cytoscape, Molecular Complex Detection (MCODE) software was used, of which the degree cut-off=2, the cut-off of node score=0.2, the k-core=2, while the max. Depth=100.

Survival analysis of hub genes

Based on 50 osteosarcoma patients in one dataset (GEO21257), in terms of patient survival, the relevance of the genes identified was displayed in the form of Kaplan-Meier plot.

Results

Identification of DEGs by meta-analysis

In accordance with the inclusion criteria, six GEO datasets as follows were obtained from the NCBI: GSE11414, GSE12865, GSE14359, GSE16102, GSE42352, GSE42572 (see the “Materials and Methods”, Figure 1A). 173 samples of osteosarcoma and 29 samples of normal tissues in total were used in meta-analysis. Using Affymetrix gene chips and Illumina gene chips, the GEO Platform (GPL) files were obtained from total six datasets (Table 1). Total 10390 common genes were identified throughout all datasets, and in accordance with the combination of p-values and Z scores, a meta-analysis of the expression profiles of multiple genes was performed on two platforms. Total 218 DEGs were identified in accordance with the combination of t-values (with the threshold being 5) and Z scores (with the limit being 4) (Figure 1B). And total 75 genes down-regulated and 142 genes upregulated were identified (see supplementary materials Table S1).

Figure 1

Screening for DEGs. (A). Selection procession of microarray datasets for meta-analysis. (B). The 218 overlapping DEGs based on |pval_test|>5 and |Z|>4 were detected using Venny2.1.0.

Table 1

Characteristic of individual studies retrieved from Gene Expression Omnibus for meta-analysis.

Dataset	Samples	Case/Control	Country	PMID	Platforms	Gene#	Gene chip
GSE11414	6	4/2	Canada	18698372	GPL6244	33297	Affymetrix Human Gene 1.0 ST Array
GSE12865	14	12/2	Canada	19286668	GPL6244	33297	Affymetrix Human Gene 1.0 ST Array
GSE14359	20	18/2	Germany	21166698	GPL96	22283	Affymetrix Human Genome U133A Array
GSE16102	51	48/3	USA	20028558	GPL96	22283	Affymetrix Human Genome U133A Array
GSE42352	99	84/15	Norway	24447333/ 23688189	GPL10295	48701	Illumina HumanWG-6 v2.0 expression beadchip
GSE42572	12	7/5	Norway	26106474	GPL13376	48701	Illumina HumanWG-6 v2.0 expression beadchip

According to the inclusion criteria, the following six GEO datasets from the NCBI were obtained: GSE11414, GSE12865, GSE14359, GSE16102, GSE42352, GSE42572.A total of 173 osteosarcoma samples and 29 normal tissue samples were used in meta-analysis. The GEO Platform Files (GPLs) from the six datasets were obtained using Affymetrix gene chips and Illumina gene chips.

Analysis of GO term and KEGG pathway enrichment

To make a deep investigation of the DEGs’ functions, DEGs were put into the category of functional GO and KEGG, and then an analysis of pathway enrichment was carried out. Consistent with the p-value (P<0.05), selecting the top ten terms among all terms enriched in the KEGG category was made. The most enriched terms regarding the DEGs in KEGG pathway were Apoptosis (hsa04210), Lysosome (hsa04142), Calcium signaling pathway (hsa04020), Cell adhesion molecules (hsa04514), Endocytosis (hsa04144) (Figure 2A). For cellular components, it’s demonstrated by the results of GO analysis that there was separate enrichment of DEGs in cytoplasmic vesicle part (GO:0044433), lysosome (GO:0005764), lytic vacuole (GO:0000323), and secretory granule (GO:0030141) (Figure 2B). Based on the above pathway analysis, we believed that osteosarcoma’s occurrence and development were closely related to macrophages. For Cellular Component analysis in GO term, the significant gene sets in the enrichment results were shown in a directed acyclic graph (DAG), and the diagram showed that the lysosomal related pathways were closely associated with the development of OS (Figure 2B). We conducted GSEA analysis based on a negative z score and get a similar conclusion (see supplementary materials Table S2-4).

Figure 2

DEGs pathway analysis showed that they were associated with macrophages. (A). KEGG enrichment analysis of DEGs indicated the lysosomal pathway. (B). DEGs of lysosomal pathway in Cellular Component analysis in GO term.

LM22 GSEA pathway analysis and EPIC analysis in the tumor microenvironment of osteosarcoma

To investigate the immune-related characteristics that might be involved in immunoreactive osteosarcoma, a GSEA enrichment prediction of the tumor microenvironment was conducted. Newman et al. constructed a leukocyte signature matrix that enabled further applications on immune infiltration by applying CIBERSORT (short for cell-type identification by estimating relative subsets of RNA transcripts), which was a new approach that allowed the RNA mixtures to be analyzed on a large scale. Based on LM22 gene sets, GSEA analysis results showed that the monocyte pathway was significantly activated (Figure 3A). Compared with the normal group, genes associated with monocytes were significantly upregulated in osteosarcoma tissues (Figure 3B). We also performed EPIC analysis to estimate the number of immune cells existing in the tumor microenvironment of osteosarcoma. EPIC was made based on a unique collection of RNA-seq reference gene expression profiles coming from the circulating immune cells or the tumor invasive non-malignant cells. We obtained four datasets in the GEO dataset, including GSE126209, GSE87624, GSE99671, and GSE57925 (Table 2). Consistent with GSE analysis, we concluded that the proportion of macrophages in osteosarcoma samples was greater than that in normal samples. The proportion of macrophages was significantly increased while the proportion of lymphocytes had no statistical changes in immunoreactive osteosarcoma compared with normal tissues.

Figure 3

Tumor microenvironment analysis revealed abnormal macrophage function in osteosarcoma. (A). the result of LM22 based GSEA pathway analysis. (B). the enrichment plot of monocytes. (C). EPIC analysis showed osteosarcoma patients had more macrophage cells than normal tissues.

Table 2

Characteristic of individual studies retrieved from Gene Expression Omnibus for EPIC.

Dataset	Samples	Case/Control	Country	PMID	Platforms	Gene#	Gene chip
GSE126209	11	6/5	China	NO	GPL20301	60493	Illumina HiSeq 4000 (Homo sapiens)
GSE87624	47	44/3	USA	29066513	GPL11154	21887	Illumina HiSeq 2000 (Homo sapiens)
GSE99671	51	33/18	Estonia	29050494/ 29250102	GPL20148	22775	AB 5500xl-W Genetic Analysis System (Homo sapiens)
GSE57925	3	3/0	USA	25961939	GPL11154	22259	Illumina HiSeq 2000 (Homo sapiens)

The following four GEO datasets from the NCBI were obtained: GSE126209, GSE87624, GSE99671, GSE57925.A total of 86 osteosarcoma samples and 26 normal tissue samples were used in EPIC analysis. The GEO Platform Files (GPLs) from the four datasets were obtained using illumine HiSeq gene chips.

Screening of hub gene and module from PPI network and the survival analysis

In the first place, DEGs’ PPI network was determined. Based on the STRING database, this PPI network was composed of 217 nodes as well as 274 edges. As for the top six hubs with degree centrality as screened from the PPI network, they were taken as the hub genes (Figure 4A). For these hub genes, they included CD209 molecule (CD209), CD34 molecule (CD34), CD86 molecule (CD86), CD163 molecule (CD163), CD80 molecule (CD80) and Fc fragment of lgG receptor IIb (FCGR2B). Most important of all, these hub genes were all upregulated in osteosarcoma tissues (Figure 4B). To make an investigation of the relationship between these genes and the patient prognosis, the hub genes were carried out a survival analysis. We found that GSE21257 contains only cases without normal. It is not used for meta-analysis, but only for survival analysis. Our study found that the expression of CD163, as the marker of macrophages, was associated with the OS of the osteosarcoma patients in a negative way (OS_HR=2.2889, P=0.0497). These results suggested that macrophage was a risk factor in patients with osteosarcoma (Figure 5A-B).

Figure 4

Establishment of PPI network. (A). Macrophage-related modules obtained from PPI network of DEGs using the MCODE software. (B). The expression of Macrophage-related modules between osteosarcoma versus normal tissue in GEO dataset (GSE42352).

Figure 5

The expression of the marker (CD163) in macrophages and its prognosis. (A). meta-analysis for comparing CD163 expression in the osteosarcoma and para-carcinoma tissues based on six GEO datasets. The forest plot shows that CD163 expression is upregulated in cancer. (B). CD163 significantly correlates with poor prognosis of OS(GEO21257).

Discussion

We collected all the osteosarcoma data sets in GEO datasets and selected the data sets meeting the requirements for meta-analysis. According to the two cut-off values of t-values and Z scores, 218 DEGS were determined. Next, it’s shown by analysis of GO term as well as KEGG pathway enrichment that macrophages played a significant role in the development of osteosarcoma. EPIC analysis indicated that macrophages were increased in osteosarcoma, and lm22 GSEA pathway analysis suggested abnormal activation of the macrophage pathway. In addition, PPI network analysis showed that the macrophage module was the core module, and survival analysis indicated that high expression of the marker of macrophage (CD163) predicted a worse prognosis in patients with osteosarcoma.

The tumor microenvironment (known as TME) is constituted by various cell types (e.g., the endothelial cells, the fibroblasts, and the immune cells) as well as extracellular components (e.g., the cytokines, the growth factors, the hormones, and the extracellular matrix) as surrounding the tumor cells and receiving nourishment from a vascular network. In terms of TME, it does not only exert an important effect on initiation, progression, and metastasis of tumors but also plays a crucial role in the therapeutic efficacy of tumors[6]. It should be noticed that the therapeutic efficacy could be supported or obstructed by the immune cells. Still, the activation status and location of the immune cells themselves in the TME could be varied. Based on EPIC, this study used the gene expression matrix of patients with osteosarcoma to find that the proportion of macrophages in the microenvironment of osteosarcoma increased compared with the normal tissues, and the survival analysis found that the change had a significant impact on the survival of the patients, which deepened the understanding of the microenvironment of osteosarcoma and enhanced the knowledge of osteosarcoma.

There has been strong evidence showing that tumor-associated macrophages (known as TAMs) form a key regulator for a therapeutic response within the TME. For bone marrow monocytes, which come from progenitor cells of bone marrow, they can enter the tumors through blood circulation and then be differentiated into macrophages later. Based on their polarization status, macrophages can be classified into subcategory M1 and subcategory M2. Subcategory M1 macrophages could be activated by Th1 cytokine interferon γ (known as IFNγ) and the microbial products. In contrast, subcategory M2 macrophages would differentiate as a consequence of Th2 cytokines, with IL-4, IL-10, IL-13, and so on included[11,12]. In the context of TAMs, it’s considered that M1 macrophages could exert an effect of tumor-destroying, while the M2 macrophages could promote the development of tumors. Being both plastic and reversible, M1 and M2 TAMs could have their functional polarization regulated by TME significantly[13,14]. Studies have shown that cd163 was the marker of M2 macrophage in bone marrow[15]. Therefore, this study first indicated that the M2 macrophage in TME of osteosarcoma increased and was harmful to osteosarcoma patients.

During the therapeutic process of cancer, a misdirected response to tissue repair with coordination by TAMs could be induced by the chemotherapeutic agents[16]. This might promote the growth of tumors and limit the anti-tumor activity. It has been demonstrated by evidence both in vitro and in vivo that the resistance to some chemotherapeutic agents (for example, the 5-fluorouracil, the doxorubicin, the gemcitabine, the paclitaxel, and the platinum compounds) and anti-VEGF treatment could be mediated by TAMs[17-21]. Myeloid-derived suppressor cells could migrate to the tumor site actively as well and then differentiate into TAMs rapidly[22]. With the enhancement of the knowledge on tumor immunity and immune escape mechanisms, some tumor immunotherapies that are new and more effective have come into existence. As for methods that destroy the immune tolerance to TME (short for tumor microenvironment), make effective prevention of tumor immune escape and stimulate antitumor immunity in the body more effectively, while not making stimulation or even suppression of the negative immune responses, they are novel design approaches for new modalities of therapy. Combined with previous studies, this study suggested that M2 macrophages cells in osteosarcoma may act as a potential new therapeutic target.

As for macrophages’ functions, they were as changeable as their lineages and had a significant effect on regulating normal homeostasis and the development of diseases[23,24]. Increasing evidence showed that macrophages played a principal role in the remodeling of normal and diseased tissues, such as angiogenesis, breakdown of basement membranes, filtration of leukocytes, suppression of immunity, and so on[23]. Therefore, with TME (tumor microenvironment), macrophage has become a drug target with a central role for various states of diseases.

As a receptor tyrosine kinase that is expressed on macrophages residing in tissues, at the sites of inflammation, RON (Recepteur d’Origine Nantais) would bind to the macrophage-stimulating protein (known as MSP) upon activation of proteolysis[25]. As for both prostate and breast cancer models, it has been shown that RON signaling in macrophages could have CD8+ T cells’ anti-tumor functions impaired, thus promoting the growth of tumors and outgrowth of metastasis, respectively[26,27]. The pharmacologic blockade against RON kinase by the inhibitor such as BMS-777607/ASLAN002 could boost the number of TNF-α-secreting macrophages pro-inflammatory and reduce outgrowth of metastasis in lungs in the breast tumor model of PyMT-MSP[27]. Therefore, drugs targeting macrophages may be a new therapeutic target for osteosarcoma.

To the best of our knowledge, this view about M2 macrophages in osteosarcoma has never been reported in the literature before, so our findings not only suggested that changes in the immune microenvironment in osteosarcoma but also may become a new direction of immunotherapy. This integrated analysis may help to identify potential targets for osteosarcoma treatment, explore relevant pathways, increase the rigor of preclinical, experimental design, and provide a new treatment method for osteosarcoma.

In conclusion, this study has systematically validated the results of the studies carried out previously and filled up the gap in OS’s field on a large-scale meta-analysis. Also, for the hub gene (CD163) and the macrophage cell capable of being used as a novel biomarker in promoting early diagnosis and development of therapeutic approaches, the evaluation was made, which provided strong evidence for making individualized treatment of OS based on genomes in the future.

Author contribution

SZ, JP and XC conceived and designed the study, and drafted the manuscript. SZ, JJ, TW and XC collected, analyzed and interpreted the data. SZ and JP revised the manuscript for important intellectual content. All authors read and approved the final manuscript.