- Open Access
The role of SPI1-TYROBP-FCER1G network in oncogenesis and prognosis of osteosarcoma, and its association with immune infiltration
BMC Cancer volume 22, Article number: 108 (2022)
Osteosarcoma is an aggressive malignant bone sarcoma worldwide. A causal gene network with specific functions underlying both the development and progression of OS was still unclear. Here we firstly identified the differentially expressed genes (DEGs) between control and OS samples, and then defined the hub genes and top clusters in the protein–protein interaction (PPI) network of these DEGs. By focusing on the hub gene TYROBP in the top 1 cluster, a conserved TYROBP co-expression network was identified. Then the effect of the network on OS overall survival was analyzed. Gene Ontology (GO), Kyoto Encyclopedia of Genes and Genomes (KEGG) enrichment analyses and Gene Set Enrichment Analysis (GSEA) were used to explore the functions of the network. XCell platform and ssGSEA algorithm were conducted to estimate the status of immune infiltration. ChEA3 platform, GSEA enrichment analysis, and Drug Pair Seeker (DPS) were used to predict the key transcription factor and its upstream signal. We identified the downregulated SPI1-TYROBP-FCER1G network in OS, which were significantly enriched in immune-related functions. We also defined a two-gene signature (SPI1/FCER1G) that can predict poorer OS overall survival and the attenuated immune infiltration when downregulated. The SPI1-TYROBP-FCER1G network were potentially initiated by transcription factor SPI1 and would lead to the upregulated CD86, MHC-II, CCL4/CXCL10/CX3CL1 and hence increased immune infiltrations. With this study, we could better explore the mechanism of OS oncogenesis and metastasis for developing new therapies.
Osteosarcoma (OS) is the most common primary malignant bone sarcoma worldwide, particularly in children and adolescents [1, 2]. OS is a highly invasive tumor with 20% of patients being diagnosed at metastatic stage and 50% of patients with lung metastases at late stage. OS was due to a dysfunctional osteoblastic differentiation from the primitive mesenchymal bone-forming cells . Now the treatment for OS has reached a plateau to the patients with metastases, the 5-year survival rate of which is under 30% . Therefore, identifying the key gene network and their molecular mechanisms underlying both the pathogenesis and progression of OS is necessary and urgent.
Some studies have been published, which focused on the role of single oncogene in OS, such as p53 (arresting cell cycle and promoting apoptosis) , SOX2 (enhancing cell stemness and migration) , MALAT1 (promoting proliferation and metastasis) , IGF-2 (affecting osteoblast differentiation) , cyclin E1 (inhibiting proliferation and increasing chemotherapeutics sensitivity) , and PTH (affecting migration) . In recent years, some studies identified a cluster of hub genes as OS biomarkers via differentially expressed genes (DEGs) analysis followed by protein–protein interaction (PPI) network analysis or weighted gene co-expression network analysis (WGCNA). In this kind of analyses, multiple clusters of hub genes have been identified as OS characteristic biomarkers or metastatic biomarkers involving functions such as focal adhesion, type I interferon signaling and cell–cell adhesion [11,12,13,14,15]. However, a conserved gene network involved in both OS oncogenesis and OS prognosis have not been fully analyzed. Furthermore, the common gene network and functions across cancer types have not been fully discussed.
Here we integrated DEG results from GSE33382, GSE12865, GSE16088 and GSE14359 datasets and analyzed the survival data from TARGET-OS and GSE21257 datasets to identify the immune-infiltration associated gene network (i.e., SPI1-TYROBP-FCER1G network) underlying both OS oncogenesis and prognosis. Firstly, we identified the DEGs and defined the hub genes and top clusters in the PPI network of those DEGs. Secondly, we did univariate Cox regression analysis to identify the 9 prognostic genes in the five top clusters. To furtherly explore the function of TYROBP-centered top 1 cluster, we did an analysis across cancer types to identify a more conserved and extensive TYROBP co-expression network, which was found to be significantly associated with leukocyte proliferation and T cell activation later. By PPI analysis and least absolute shrinkage and selection operator (LASSO) Cox regression analysis on the TYROBP co-expression network, the two-gene signature of FCER1G and SPI1 were identified as a novel prognosis biomarker. Furthermore, the transcriptional levels of FCER1G and SPI1 can also distinguish the OS patients into two clusters with different prognosis and immune cell infiltrations, suggesting a potential treatment target for OS and other TYROBP-associated cancers.
Materials and methods
Data resource and processing
Four microarray datasets of OS samples were collected from GEO database, including GSE33382, GSE12865, GSE16088 and GSE14359. The selection criteria included: 1. including both normal control and osteosarcoma tissues; 2. from Homo sapiens; 3. more than 10 samples. The platform annotation documents were also downloaded from GEO and annotated for microarray probes by “merge” R command.
Both expressional data and survival data from TARGET-OS and GSE21257 datasets were downloaded for analyzing the effect of gene expression on prognosis (1. including both osteosarcoma tissue and survival data; 2. from Homo sapiens; 3. more than 20 samples).
TCGA Pan-Cancer TPM data (TOIL workflow processed ) were obtained from UCSC Xena Browser (https://xenabrowser.net/datapages/). The correlations between the expressional levels of TYROBP and other genes were assessed using Spearman correlation analysis (R > 0.6 and p < 0.05 as significance). GEPIA2 platform was used to analyze the correlation between SPI1, TYROBP and FCER1G across multiple cancer types .
ENCODE PU.1 ChIP-seq peaks data from multiple cell lines (i.e., GM12878, GM12891, K562, HL-60) were downloaded  and visualized by software IGV . RNA-seq and ChIP-seq data from THP-1 cell line were downloaded from GSE69284 and GSE89178, the peaks of which were called by HOMER . ChIP-X Enrichment Analysis 3 (ChEA3) platform was used for transcription factor (TF) prediction by transcription factor enrichment analysis that ranks TFs associated with user-submitted gene sets .
Analysis of differentially expressed genes
GEO2R online tool was used to systematically measure the DEGs in GSE33382, GSE12865 and GSE14359 (between control and OS). ‘GEOquery’ R package and ‘limma’ R package were used to measure the DEGs in GSE16088 (between control and OS) [22, 23]. Significant DEGs were selected by Benjamini & Hochberg (BH, False discovery rate) adjusted p < 0.05 and |log2 fold change (FC)|> 0.58. Then the overlapped results were visualized by Venn diagram via “ggplot2” R package .
Functional enrichment analysis
Gene Ontology Analysis on Gene Clusters Gene ontology (GO) (including biological process (BP), molecular functions (MF), and cellular components (CC) terms)  and Kyoto Encyclopedia of Genes and Genomes (KEGG) pathways  are used to annotate the functions of genes via “clusterProfiler” R package  and virtualized by “ggplot2” R package.
Gene set enrichment analysis (GSEA) was performed by pre-ranking genes based on R value of TYROBP correlation level. We subsequently run the GSEA preranked analysis with the c2.cp.reactome.v7.2.symbols.gmt (reactome), h.all.v7.2.symbols.gmt (signaling pathway), and c3.all.v7.2.symbols.gmt (transcriptional factor motif) gene set from MsigDB  via “clusterProfiler” R package.
Protein–protein interaction network
The Search Tool for the Retrieval of Interacting Genes (STRING) database is an online tool that is used to develop protein–protein interaction (PPI) networks . The network data of proteins from DEGs were downloaded from STRING to visualize the protein interaction network in the Cytoscape software . In the network of DEGs, the 5 top clusters (highly interconnected regions) were identified by MCODE plugin in Cytoscape, which can find clusters (i.e., the highly interconnected regions) by vertex weighting and isolating the dense regions in a network . The degree of gene connection was analyzed by Cytoscape and reflected by its area. In the TYROBP network, we identified the ranked hub genes by cytoHubba plugin (with Maximal Clique Centrality (MCC) algorithm)  in Cytoscape.
Survival analyses were performed using univariate Cox regression and LASSO regression via “survival”  and “glmnet”  R package. For LASSO regression, the final formula of the risk signature was given as risk score = expression level of Gene1*β1 + expression level of Gene2*β2 + … + expression level of Genen* βn, in which β represents the regression coefficient of each variable. The receiver operating characteristic (ROC) curves for the imaging markers were constructed for overall survival, and the areas under the ROC curves (AUC) were estimated empirically with the trapezoid rule with “pROC”  and “ggplot2” R package.
Immune infiltration analysis
Platform xCell (xCell signature)  and single-sample gene set enrichment analysis (ssGSEA, based on immune signature by “gsva” R package ) were used to reckon the proportion and score of immune cell infiltration in OS tumor environment, respectively. We filtered out the cell types with proportion less than 0.01 in xCell results. Boxplots of proportions and normalized scores were performed with “ggplot2” R package. Wilcox test was performed for the comparison between groups. ssGSEA, as an extension of GSEA, can calculate separate enrichment scores of immune cells for the pairing of OS samples and the gene set of immune cell signatures. XCell can calculate separate proportions of immune cells in microenvironment based on gene signatures for 64 cell types.
The DEGs from multiple OS datasets and the immune-related function of downregulated genes
Four of GEO datasets with both normal control and osteosarcoma tissues were included in the study. It’s notable that there were a few of differences between these datasets, reflecting the complementarity between them. Firstly, in GSE12865 and GSE14359, the control samples were normal primary osteoblasts, while in GSE33382 the control samples were osteoblasts derived from osteogenic differentiated bone-marrow-derived mesenchymal stem cells. In GSE16088, the control samples were from an unspecific osteoblast cell line H-012706. Furthermore, although OS samples were all OS tissues not cell line in the four datasets, in GSE33382 the OS samples were from high-grade OS pre-chemotherapy biopsy. In GSE12865, OS samples were all from pediatric OS patients. In GSE14359, OS group included a part of lung metastasis OS tissues. For a comprehensive analysis, the DEGs ( |log2(Fold Change)|> 0.58 and adjusted p < 0.05) were identified in 4 OS datasets independently (Supplementary Fig. 1A). The overlapped DEGs in > = 3 datasets were considered as credible DEGs for the next functional enrichment analysis (Fig. 1A-B). In the end, we found 124 of upregulated genes (Fig. 1A) and 98 of downregulated genes (Fig. 1B).
Among them, the 124 upregulated genes were enriched in annotations such as collagen-containing extracellular matrix, smooth muscle cell proliferation, cell-substrate adhesion and p53 signaling pathway (Fig. 1C), while the 98 downregulated genes were enriched in immune-related functions including MHC class II protein complex/complex binding, lymphocyte/leukocyte proliferation, and so on, indicating a significant loss of MHC class II protein complex and immune cell activation ability in OS samples compared with that in control samples (Fig. 1D).
The hub genes and top clusters in PPI network of DEGs, and the immune-related function of top 1 cluster genes
To define the hub genes and functional gene clusters in DEGs, we put upregulated and downregulated genes together to construct the PPI network by Cytoscape based on STRING database. According to the connection degrees of DEGs, we found the hub genes with the top 10 degrees (reflected by their areas in Fig. 2A) in the network (Supplementary Fig. 1B) including FN1, TYROBP, EGFR, CSF1R,C1QB, PLEK, LCP2, C1QA, CD86 and VWF, which were enriched in functions including osteoclast differentiation/integrin binding/collagen-containing extracellular matrix/glial cell activation (Fig. 2B), and potentially initiated the whole PPI network as the hub genes.
Furthermore, we defined the top 5 clusters by MCODE plugin in Cytoscape and analyzed their functions respectively (Fig. 2C-G). Notably, genes from top 1 cluster were all downregulated in OS and enriched in immune-related functions such as positive regulation of immune cell activation/cytokine production, and MHCII class II protein complex (Fig. 2C). Top 2 cluster genes were about the functions of lipoprotein particle binding, granule membrane and phagocytosis (Fig. 2D). Top 3 cluster genes were enriched in functions of actin cytoskeleton reorganization (Fig. 2E). Top 4 cluster genes were mainly associated with fibroblast proliferation (Fig. 2F). Top 5 cluster genes were most significantly associated with cell migration including smooth muscle cell migration and leukocyte migration (Fig. 2G).
Five genes from top 1 cluster can also predict patient prognosis
To furtherly explore the potential effect of the 5 top clusters on OS patient prognosis, we did univariate Cox regression based on TARGET-OS database and found 9 of the top cluster genes can also predict the prognosis of OS patients (Fig. 3A-B). Interestingly, 5 of the 9 prognostic genes, i.e., FCER1G, CSF1R, C1QA, TYROBP, and C1QB, were from the TYROBP-centered top 1 cluster, suggesting the importance of top 1 cluster in OS progression (Fig. 3B). All the 5 genes from top 1 cluster, including TYROBP itself (Fig. 3C), had a Hazard Ratio (HR) less than 1 (Fig. 3B), suggesting their protective effect on good OS prognosis (i.e., the higher expression associated with good prognosis). By GSEA enrichment analysis, the TYROBP co-expressed genes (|R|> 0. and p < 0.05) in OS were enriched in annotations including innate/adaptive immune system, cytokine signaling, neutrophil degranulation, and TYROBP causal network (Fig. 3D), indicating the strong connection between TYROBP and immune functions. Moreover, although FN1 had the highest degree in the whole network, it did not significantly impact the OS prognosis (Fig. 3E).
The conserved TYROBP co-expression network across multiple cancer types and its classic immune activation function
As the second highest degree gene in DEGs network and the highest degree gene in top 1 cluster, TYROBP also significantly impact the prognosis of OS patients (Fig. 3C). Therefore, we furtherly identified the conserved TYROBP co-expressed network across multiple cancer types to explore its classic function in tumor progression. Firstly, we found that the expressional level of TYROBP was significantly (adjusted p < 0.05) lower in cancers including READ, PAAD, LUAD, LUSC, COAD, and OS, while the expressional level was higher in cancers including THCA, BRCA, STAD, ESCA, and GBM (Fig. 4A). Then we integrated all the down-regulated co-expression genes of TYROBP in READ, PAAD, LUAD, LUSC, COAD, and OS, and all the up-regulated co-expression genes in THCA, BRCA, STAD, ESCA, and GBM (Fig. 4B). By overlapping the 80 of down-regulated co-expression genes and 109 up-regulated co-expression genes, we got 54 of conversed TYROBP co-expressed genes (R > 0.6 and p < 0.05) across multiple cancer types (Fig. 4B). These genes were enriched in functions such as mononuclear cell/lymphocyte/leukocyte proliferation, T cell activation, and antigen processing and presentation (Fig. 4C), suggesting a classic immune activation function and a potential immune infiltration-associated function.
A two-gene signature (FCER1G and SPI1) from TYROBP co-expression network can predict OS prognosis
To explore the effect of the conserved TYROBP co-expression network (Supplementary Fig. 1C) on OS prognosis, we identified the top 10 hub genes in the network by PPI analysis and defined a two-gene signature from the 10 hub genes by the overall survival-based LASSO Cox regression model. Firstly, we found the 10 hub genes including SPI1 (first rank), TYROBP (second rank), FCER1G (third rank), ITGB2, C1QB, C1QA, LY86, LCP2, CCR1, and AIF1 (Fig. 5A). Even across the ten cancer types in Fig. 4B, the expressional levels of SPI1, TYROBP, FCER1G were also highly correlated with R > 0.8 and p = 0 (Fig. 5B). Therefore, we also name the TYROBP co-expression network as SPI1-TYROBP-FCER1G network. Furthermore, a LASSO Cox regression analysis identified an optimal two-gene predictive signature (Fig. 5C) with risk score as follows: (-0.1405) × FCER1G expression value + (-0.0282) × SPI1 expression value. The AUC of risk score was 0.673 as independent prognostic indicator (Fig. 5D), which was higher than the AUCs (Fig. 5E) of FCER1G (0.664), TYROBP (0.634), and SPI1 (0.658) expression value respectively. The overall survival of OS patients with higher risk score is significantly shorter than those with lower risk score (cut off by median risk score) from overall survival curve (Fig. 5F).
Different levels of immune infiltration between the two clusters of OS patients based on the expressional levels of FCER1G and SPI1
Because dividing OS patients by median risk score is relative arbitrary, we divided OS patients according to the expressional levels of FCER1G and SPI1 in TARGET-OS database by the consensus clustering analysis. We found cluster 1, which had lower expressional levels of FCER1G and SPI1, was significantly associated with poorer overall survival (p = 0.010, Fig. 6A). In addition, cluster 1 also had lower stromal and immune scores (p < 0.01 and p < 0,001 separately, Fig. 6B), compared with cluster 2. To explore the alteration of different immune cell subsets, we estimate the immune cell infiltrations via both ssGSEA algorithm and xCell platform. We found almost all immune cells had higher scores in cluster 2 than that in cluster 1 in addition to CD56dim NK cells and Th2 cells (Fig. 6C). Besides, there were higher proportions of CD4 + Tem cells, endothelial cells (including ly and mv endothelial cells), macrophages (including M1 and M2), monocytes, dendritic cells (including aDC, cDC, iDC and pDC) in cluster 2 than that in cluster 1 (Fig. 6D), after filtering cell subset less than 1% (0.01). Furthermore, xCell also identified the significantly higher immune score and microenvironment score in cluster 2 than that in cluster 1, suggesting more immune infiltration and lower tumor microenvironment abundance (Fig. 6E).
To validate the result, we did the same analyses on GSE21257 data. Similarly, cluster 1 had significantly poorer overall survival compared with cluster 2 (p = 0.021, Fig. 6F-G). Almost all immune cells had higher scores in cluster 2 than that in cluster 1 in addition to CD56dim NK cells and pDC cells (Fig. 6H). Besides, there were higher proportions of CD4+ T cells (including T cells and Tem), CD8+ T cells (including T cells, Tcm and Tem), B cells (including B cells and class-switched memory B cells), endothelial cells, macrophages (including M1 and M2), monocytes, dendritic cells (including iDC and pDC), basophils, eosinophils, erythrocytes, DC, and iDC in cluster 2 than that in cluster 1 (Fig. 6I), after filtering cell subset less than 1%. However, there were lower proportions of plasma cells and smooth muscle cells in cluster 2 than cluster 1 (Fig. 6I). Similarly, xCell also identified the significantly higher immune score/microenvironment score but lower stroma score in cluster 2 than that in cluster 1, indicating higher immune infiltration but lower tumor purity and stroma cell proportion (Fig. 6J).
These results indicated that low expressional levels of FCER1G and SPI1 could predict an attenuated immune infiltration and poorer prognosis, suggesting a potential immune escape and metastasis tendency. By comparing the expressional level between OS patients with metastases after diagnosis of primary tumor and patients without metastases in GSE21257, we found the expressional levels of TYROBP (p = 0.0003, Log2FC = -1.16), FCER1R (p = 0.0002, Log2FC = -1.19) and SPI1 (p = 0.0059, Log2FC = -0.599) were all down-regulated in metastases group, consistent with the above result.
Transcription factor SPI1 may initiate the expression of whole TYROBP network and be impacted by upstream TNF-α
Although TYROBP is in the center of its whole co-expression network, it cannot regulate gene expression as a transcription factor. Therefore, we predicted the key transcription factor of the 55 genes of the TYROBP network by ChEA3 platform. SPI1 ranked second in the top 50 predicted transcription factors. Furthermore, in addition to the 5 genes (i.e., C1orf162, SLAMF8, HLA-DPB1, HLA-DRA, and HLA-DRB1), all other 50 genes in the network, including TYROBP, FCER1G, and SPI1 itself, can be bound and regulated by SPI1 (Fig. 7A). Interestingly, GSEA enrichment analysis also showed that transcription factor PU.1/SPI1 motifs were enriched in TYROBP positively correlated genes (Fig. 7B). ChIP-seq and RNA-seq data validated the binding of SPI1 on TYROBP (Fig. 7C) and other 9 of network hub genes (data not shown). Notably, the low expressional level of SPI1 was significantly associated with poorer overall survival in TARGET-OS dataset (Fig. 7D).
To identifying the upstream factor of SPI1 in OS, by GSEA analysis, we found ‘TNF-α signaling via NF-κB’, ‘IFN-α/γ response’, ‘inflammatory response’ and ‘IL2-STAT5 signaling’ pathways were all enriched in TYROBP positively correlated genes (Fig. 7E). Moreover, by the analysis of Drug Pair Seeker (DPS), we found TNF-α and IGF2 can partly rescue the under-expression of TYROBP network and the over-expression of TYROBP negatively-correlated genes in OS (Fig. 7F). These results suggested the importance of TNF-α in the activation of SPI1 and TYROBP network, and a potential deficiency of tumor microenvironmental TNF-α in OS and OS with poorer prognosis.
The TYROBP network-associated immune molecules affecting immune infiltration and activation
To furtherly explore the mechanism of downregulated immune infiltrations in OS, we overlapped the immune stimulators, MHC molecules, chemokines, and immune inhibitors with 965 of TYROBP positively-correlated genes (R > 0.4, p < 0.05) in TARGET-OS dataset and 491 of down-regulated genes in > = 2 GEO datasets (i.e., GSE33382, GSE12865, GSE16088 and GSE14359). We identified 1 of immune stimulators (CD86), 5 of MHC molecules (HLA-DMA, HLA-DMB, HLA-DPA1, HLA-DPB1, HLA-DRA), 3 of chemokines (CCL4, CXCL10, CX3CL1), and 2 of immune inhibitors (CSF1R, HAVCR2), which were both TYROBP positively-correlated and down-regulated in OS samples compared with that in normal samples (Fig. 8A-B), which potentially underlay the mechanism of attenuated immune infiltrations by the under-expression of SPI1-TYROBP-FCER1G network in OS (Fig. 8C).
This study identified a conserved TYROBP co-expression network that was associated with both OS oncogenesis and prognosis. The whole workflow was summarized in Fig. 8D. Firstly, we found five functional gene clusters of DEGs between control and OS samples. Then we focused on the TYROBP-centered top 1 cluster that was down-regulated in OS and associated with MHC II complex, immune cell activation and cytokine production. The under-expression of five genes (i.e., FCER1G, CSFR1, C1QA, TYROBP, and C1QB) from top 1 cluster were also associated with poorer overall survival of OS patients. Furthermore, we extended the top 1 cluster into a more conserved TYROBP co-expression network with 10 hub genes (SPI1, TYROBP, FCER1G, ITGB2, C1QB, C1QA, LY86, LCP2, CCR1, AIF1) by an analysis across cancer types, which was associated with immune cell proliferation and T cell activation. Then, the risk score combining the expression levels of both SPI1 and FCER1G were constructed to predict the overall survival of OS patients. Moreover, we found the distinguished immune infiltrations between two clusters of OS patients with different overall survival times, suggesting the importance of SPI1-TYROBP-FCER1G network in upregulating the microenvironmental immune infiltration and promoting better prognosis. In addition, we identified SPI1 as the key transcription factor in the network, and found it’s potentially associated with the upstream TNF-α signaling. We also identified the TYROBP positively-correlated and downregulated immune molecules in OS, suggesting the immune-related mechanism of attenuated immune infiltration and activation in OS.
As a conserved and immune-related network across multiple cancer types, some molecules in the SPI1-TYROBP(i.e., KARAP or DAP12)-FCER1G network have been identified as prognostic biomarkers before in other cancer types [38,39,40,41] and OS [42,43,44]. Although the network was upregulated in certain cancer types (e.g., clear cell renal cell carcinoma, breast cancer, gastric cancer, and low-grade glioma) and positively associated with the poorer prognosis and higher immune infiltrations, the network was downregulated in OS, positively associated with higher immune infiltrations in OS, and negatively associated with the poorer prognosis of OS. Therefore, despite the same positive-correlation between TYROBP network and higher immune infiltrations across multiple cancer types [38,39,40,41,42,43,44,45], the associations between TYROBP network and overall survival in different cancer types were different, suggesting the altered role of immune infiltration or altered major infiltrating cell types among different cancer types. In a word, SPI1-TYROBP-FCER1G network is a conserved immune-related network underlying both the oncogenesis and the prognosis of OS and other cancer types.
Across multiple cancer studies, the TYROBP network was always positively-associated with immune infiltration, such as the increased infiltration of CD8+ T cells, DCs, and macrophages [38,39,40,41,42,43,44,45]. According to our study on its potential mechanism, the overexpression of TYROBP network may activate the immune infiltration by the co-expressed immune stimulators, MHC molecules, and chemokines. After recruiting immune cells into the tumor microenvironment by secreting CCL4 (collaborating with CCR5 to recruit CD8 + T cells, NK cells, monocytes, etc. [46, 47]), CXCL10 (collaborating with CXCR3 to recruit CD8 + T cells, NK cells, etc. [48,49,50]), and CX3CL1 (not only collaborating with CX3CR1 to recruit CD4 + T cells, CD8 + T cells, dendritic cells [51, 52], but also inducing OS metastasis via ICAM-1 ), molecules CD86 (co-stimulatory molecule CD86 polymorphism is also associated with increased susceptibility to osteosarcoma ) and MHCs (HLA-DPA1/HLA-DPB1/HLA-DRA) on the surfaces of the tumor cells, endothelial cells or DCs/macrophages can furtherly activated T cells in the microenvironment. Therefore, the under-expression of TYROBP network can lead to the attenuated immune infiltrations of OS due to the lack of chemokines, stimulators and MHCs, which can cause the immune escape and then the metastasis of OS. Furthermore, we also found the down-regulation of immune checkpoint inhibitors CSF1R and HAVCR2 in OS. It has been found that OS tumor cell-intrinsic CSF1R can enhances cell proliferation and metastasis by CSF-1R/JAG1 axis or ERK signaling pathway [54, 55]. HAVCR2 (i.e., TIM-3), as an immune checkpoint inhibitor, can also promote the OS progression by inducing T cell exhaustion and impair T cell function when expressed on T cells [56, 57] and inducing M2 macrophage polarization when expressed on macrophages [58, 59]. The relationship between the under-expression of CSF1R and immune infiltration in OS needs further study in the future.
As transmembrane immune signaling adaptor proteins, TYROBP and FCER1G may mediate an intracellular signal via immunoreceptor tyrosine-based activation motif (ITAM) to up-regulate the expression of SPI1 and then initiate the whole SPI1-TYROBP-FCER1G network for a competent immunological surveillance function with upregulated immune infiltrations in tumor microenvironment [60,61,62,63]. SPI1, as an oncogene, is also an important transcription factor for monocyte and macrophage identity. However, it can also be regulated by multiple signals, such as NF-κB signaling from upstream TYROBP, FCER1G [62, 63], TNF-α [64, 65] or IFN-γ [66, 67]. SPI1, as the key transcription factor for the network, may shape a positive-feedback network with the signaling from TYROBP/FCER1G or environmental TNF-α/IFN-γ, then furtherly initiate the whole SPI1-TYROBP-FCER1G network to regulate the immune infiltration and cancer prognosis (Fig. 8C). In OS, the under-expression of SPI1-TYROBP-FCER1G network in tumor cells or microenvironmental endothelial cells/myeloid cells can potentially attenuate the immune infiltration and T cell activation, leading to the metastasis/poorer prognosis of OS (Fig. 8C).
This study has several limitations. First, the relatively small number of DEGs that appeared in > = 3 datasets may be due to the heterogeneity between these GEO datasets that we have mentioned in the result, which can reduce the efficiency of finding more functional gene clusters. However, it did not affect the validity of already identified results. In addition, although current studies have obtained useful findings, more laboratory studies both in vitro and in vivo are needed to validate these bioinformatic results in the future, such as the levels of SPI1-TYROBP-FCER1G network-regulated immune molecules in OS tumor cells. Thirdly, in addition to OS tumor cells, the expressional levels of SPI1-TYROBP-FCER1G network in other potential cell types, such as myeloid cells and endothelial cells, in the microenvironment are still unclear. Furthermore, the reason why SPI1-TYROBP-FCER1G network was upregulated in certain cancer types but downregulated in the others including OS, and its potential immunological mechanism, cannot be answered by our study and need further explorations in the future.
In a word, our study identified a conserved SPI1-TYROBP-FCER1G network across cancer types for the first time, and explored its immune-related mechanism in OS, which could be used as a treatment target in OS with TNF-α or TYROBP/FCER1G stimulators to rescue the under-expression of SPI1-TYROBP-FCER1G network in OS. The drug seeker analysis also supported a standpoint  that anti-TNF-α treatment in OS may cause metastasis and poorer prognosis by sustaining the under-expression of SPI1-TYROBP-FCER1G network. We provided an essential bioinformatic mechanism foundation and a novel therapeutic target of OS oncogenesis and prognosis for future explorations.
Availability of data and materials
Publicly available datasets were analyzed in this study. The data can be found here: https://www.ncbi.nlm.nih.gov/geo/query/acc.cgi (GEO datasets), https://xenabrowser.net/datapages/ (TCGA datasets), and https://ocg.cancer.gov/programs/target/data-matrix (TARGET-OS dataset).
Activated dendritic cells
Conventional dendritic cells
Immature dendritic cells
Plasmacytoid dendritic cells
Central memory T-cells
Effector memory T cells
T cells gamma delta
Common myeloid progenitors
Mesenchymal stem cells
Common lymphoid progenitors
- ly endothelial:
Lymphatic endothelial cells
- mv endothelial:
Microvascular endothelial cells
Breast invasive carcinoma
Lung squamous cell carcinoma
Raymond AK, Jaffe N. Osteosarcoma multidisciplinary approach to the management from the pathologist’s perspective. Cancer Treat Res. 2009;152:63–84.
Zhao X, et al. Osteosarcoma: a review of current and future therapeutic approaches. Biomed Eng Online. 2021;20(1):24.
Tang N, et al. Osteosarcoma Development and Stem Cell Differentiation. Clin Orthop Relat Res. 2008;466:2114–30.
Ferguson JL, Turner SP. Bone Cancer: Diagnosis and Treatment Principles. Am Fam Physician. 2018;98(4):205–13.
Chen Z, et al. TP53 Mutations and Survival in Osteosarcoma Patients: A Meta-Analysis of Published Data. Dis Markers. 2016;2016:4639575.
Maurizi G, et al. Sox2 is required for tumor development and cancer cell proliferation in osteosarcoma. Oncogene. 2018;37:4626–32.
Duan G, et al. Knockdown of MALAT1 inhibits osteosarcoma progression via regulating the miR-34a/cyclin D1 axis. Int J Oncol. 2019;54(1):17–28.
Shimizu T, et al. IGF2 preserves osteosarcoma cell survival by creating an autophagic state of dormancy that protects cells against chemotherapeutic stress. Cancer Res. 2014;74(22):6531–41.
Wei R, et al. Cyclin E1 is a prognostic biomarker and potential therapeutic target in osteosarcoma. J Orthop Res. 2020;38(9):1952–64.
Nikitovic D, et al. Parathyroid hormone/parathyroid hormone-related peptide regulate osteosarcoma cell functions: Focus on the extracellular matrix (Review). Oncol Rep. 2016;36(4):1787–92.
Jin Z, et al. Cross-Species Gene Expression Analysis Reveals Gene Modules Implicated in Human Osteosarcoma. Front Genet. 2019;10:697.
Niu J, et al. Identification of Potential Therapeutic Targets and Immune Cell Infiltration Characteristics in Osteosarcoma Using Bioinformatics Strategy. Front Oncol. 2020;10:1628.
Yao N, et al. Identification of potential crucial genes associated with vasculogenic mimicry in human osteosarcoma based on gene expression profile. Neoplasma. 2020;67(2):286–95.
Jia Y, et al. Identification of potential gene signatures associated with osteosarcoma by integrated bioinformatics analysis. PeerJ. 2021;9:e11496.
Liu J, et al. Identification of potential crucial genes and key pathways in osteosarcoma. Hereditas. 2020;157(1):29.
Vivian J, et al. Toil enables reproducible, open source, big biomedical data analyses. Nat Biotechnol. 2017;35(4):314–6.
Tang Z, et al. GEPIA2: an enhanced web server for large-scale expression profiling and interactive analysis. Nucleic Acids Res. 2019;47(W1):W556-w560.
Davis CA, et al. The Encyclopedia of DNA elements (ENCODE): data portal update. Nucleic Acids Res. 2018;46(D1):D794-d801.
Robinson JT, et al. Integrative genomics viewer. Nat Biotechnol. 2011;29(1):24–6.
Heinz S, et al. Simple combinations of lineage-determining transcription factors prime cis-regulatory elements required for macrophage and B cell identities. Mol Cell. 2010;38(4):576–89.
Keenan AB, et al. ChEA3: transcription factor enrichment analysis by orthogonal omics integration. Nucleic Acids Res. 2019;47(W1):W212–24.
Davis S, Meltzer PS. GEOquery: a bridge between the Gene Expression Omnibus (GEO) and BioConductor. Bioinformatics. 2007;23(14):1846–7.
Ritchie ME, et al. limma powers differential expression analyses for RNA-sequencing and microarray studies. Nucleic Acids Res. 2015;43(7):e47–e47.
Gómez-Rubio V. ggplot2-elegant graphics for data analysis. J Stat Softw. 2017;77:1–3.
Ashburner M, et al. Gene ontology: tool for the unification of biology The Gene Ontology Consortium. Nat Genet. 2000;25(1):25–9.
Kanehisa M, Goto S. KEGG: kyoto encyclopedia of genes and genomes. Nucleic Acids Res. 2000;28(1):27–30.
Yu G, et al. clusterProfiler: an R package for comparing biological themes among gene clusters. Omics. 2012;16(5):284–7.
Subramanian A, et al. Gene set enrichment analysis: A knowledge-based approach for interpreting genome-wide expression profiles. Proc Natl Acad Sci. 2005;102(43):15545–50.
Szklarczyk D, et al. STRING v11: protein–protein association networks with increased coverage, supporting functional discovery in genome-wide experimental datasets. Nucleic Acids Res. 2019;47(D1):D607–13.
Shannon P, et al. Cytoscape: a software environment for integrated models of biomolecular interaction networks. Genome Res. 2003;13(11):2498–504.
Bader GD, Hogue CWV. An automated method for finding molecular complexes in large protein interaction networks. BMC Bioinformatics. 2003;4(1):2.
Chin C-H, et al. cytoHubba: identifying hub objects and sub-networks from complex interactome. BMC Syst Biol. 2014;8 Suppl 4(Suppl 4):S11–S11.
Therneau TM, Lumley T. Package ‘survival.’ R Top Doc. 2015;128(10):28–33.
Simon N, et al. Regularization Paths for Cox’s Proportional Hazards Model via Coordinate Descent. J Stat Softw. 2011;39(5):1–13.
Robin X, et al. pROC: An open-source package for R and S+ to analyze and compare ROC curves. BMC Bioinformatics. 2011;12:77.
Aran D, Hu Z, Butte AJ. xCell: digitally portraying the tissue cellular heterogeneity landscape. Genome Biol. 2017;18(1):220.
Hänzelmann S, Castelo R, Guinney J. GSVA: gene set variation analysis for microarray and RNA-seq data. BMC Bioinformatics. 2013;14:7.
Jiang J, et al. Identification of TYROBP and C1QB as Two Novel Key Genes With Prognostic Value in Gastric Cancer by Network Analysis. Front Oncol. 2020;10:1765.
Lu J, et al. Elevated TYROBP expression predicts poor prognosis and high tumor immune infiltration in patients with low-grade glioma. BMC Cancer. 2021;21(1):723.
Shabo I, et al. Breast cancer expression of DAP12 is associated with skeletal and liver metastases and poor survival. Clin Breast Cancer. 2013;13(5):371–7.
Wu P, et al. TYROBP is a potential prognostic biomarker of clear cell renal cell carcinoma. FEBS Open Bio. 2020;10(12):2588–604.
Huang H, et al. Prognostic Implications of the Complement Protein C1Q and Its Correlation with Immune Infiltrates in Osteosarcoma. Onco Targets Ther. 2021;14:1737–51.
Li W, et al. Ten-gene signature reveals the significance of clinical prognosis and immuno-correlation of osteosarcoma and study on novel skeleton inhibitors regarding MMP9. Cancer Cell Int. 2021;21(1):377.
Liang T, et al. TYROBP, TLR4 and ITGAM regulated macrophages polarization and immune checkpoints expression in osteosarcoma. Sci Rep. 2021;11(1):19315.
Chen H, et al. Three-gene prognostic biomarkers for seminoma identified by weighted gene co-expression network analysis. PLoS One. 2020;15(10):e0240943.
Mukaida N, Sasaki SI, Baba T. CCL4 Signaling in the Tumor Microenvironment. Adv Exp Med Biol. 2020;1231:23–32.
Yang H, et al. A comprehensive analysis of immune infiltration in the tumor microenvironment of osteosarcoma. Cancer Med. 2021;10(16):5696–711.
Flores RJ, et al. A novel prognostic model for osteosarcoma using circulating CXCL10 and FLT3LG. Cancer. 2017;123(1):144–54.
Pradelli E, et al. Antagonism of chemokine receptor CXCR3 inhibits osteosarcoma metastasis to lungs. Int J Cancer. 2009;125(11):2586–94.
House IG, et al. Macrophage-Derived CXCL9 and CXCL10 Are Required for Antitumor Immune Responses Following Immune Checkpoint Blockade. Clin Cancer Res. 2020;26(2):487–504.
Conroy MJ, Lysaght J. CX3CL1 signaling in the tumor microenvironment. Tumor Microenvironment. 2020;1231:1–12.
Rivas-Fuentes S, et al. Regulation and biological functions of the CX3CL1-CX3CR1 axis and its relevance in solid cancer: A mini-review. J Cancer. 2021;12(2):571–83.
Wang W, et al. CD86 + 1057G/A polymorphism and susceptibility to osteosarcoma. DNA Cell Biol. 2011;30(11):925–9.
Wen ZQ, et al. Osteosarcoma cell-intrinsic colony stimulating factor-1 receptor functions to promote tumor cell metastasis through JAG1 signaling. Am J Cancer Res. 2017;7(4):801–15.
Smeester BA, et al. PLX3397 treatment inhibits constitutive CSF1R-induced oncogenic ERK signaling, reduces tumor growth, and metastatic burden in osteosarcoma. Bone. 2020;136:115353.
Anderson AC. Tim-3, a negative regulator of anti-tumor immunity. Curr Opin Immunol. 2012;24(2):213–6.
Ge W, et al. Tim-3 as a diagnostic and prognostic biomarker of osteosarcoma. Tumour Biol. 2017;39(7):1010428317715643.
Cheng Z, et al. Tumor-derived Exosomes Induced M2 Macrophage Polarization and Promoted the Metastasis of Osteosarcoma Cells Through Tim-3. Arch Med Res. 2021;52(2):200–10.
Holderried TA, et al. Molecular and immune correlates of TIM-3 (HAVCR2) and galectin 9 (LGALS9) mRNA expression and DNA methylation in melanoma. Clin Epigenetics. 2019;11(1):1–15.
Dang D, et al. Computational Approach to Identifying Universal Macrophage Biomarkers. Front Physiol. 2020;11:275.
Pan Y-G, et al. FcεRI γ-Chain Negatively Modulates Dectin-1 Responses in Dendritic Cells. Frontiers in Immunology. 2017;8:1424.
Hamerman JA, et al. The expanding roles of ITAM adapters FcRgamma and DAP12 in myeloid cells. Immunol Rev. 2009;232(1):42–58.
Wang L, et al. Identification of TYROBP and FCER1G as Key Genes with Prognostic Value in Clear Cell Renal Cell Carcinoma by Bioinformatics Analysis. Biochem Genet. 2021;59(5):1278–94.
Nutt SL. Directing the conductor: TNF regulation of HSCs. Blood, The Journal of the American Society of Hematology. 2019;133(8):771–3.
Grigorakaki C, et al. Tumor necrosis factor alpha-mediated inhibition of erythropoiesis involves GATA-1/GATA-2 balance impairment and PU.1 over-expression. Biochem Pharmacol. 2011;82(2):156–66.
Shackelford R, Adams DO, Johnson SP. IFN-gamma and lipopolysaccharide induce DNA binding of transcription factor PU.1 in murine tissue macrophages. J Immunol. 1995;154(3):1374–82.
Libregts SF, et al. Chronic IFN-γ production in mice induces anemia by reducing erythrocyte life span and inhibiting erythropoiesis through an IRF-1/PU. 1 axis. Blood. 2011;118(9):2578–88.
Kirkham B, et al. Tumor necrosis factor-alpha inhibitors: An overview of adverse effects. Waltham, MA: UpToDate; 2016.
This work was supported by the grant from the National Natural Science Foundation of China (81901656 to ML2), the Natural Science Foundation of Anhui Province (2008085QH370 to ML2), and Scientific Research of BSKY from Anhui Medical University (XJ201935 to ML2).
Ethics approval and consent to participate
All methods were carried out in accordance with relevant guidelines and regulations (Declaration of Helsinki).
Consent for publication
The authors declare that they have no competing interests.
Springer Nature remains neutral with regard to jurisdictional claims in published maps and institutional affiliations.
Differentially expressed genes, DEGs network, and SPI1-TYROBP-FCER1G network. A. The volcanos of DEGs in four OS datasets. The X-axis represented the value of log2 (fold change), while the Y-axis represented the value of -log10 (adjusted p value); the red node represented upregulated genes, while the blue node represented the downregulated genes. B. The PPI network of DEGs. The area represented the degree of each node; the red node represented upregulated genes, while the blue node represented the downregulated genes. C. The SPI1-TYROBP-FCER1G network. The area represented the degree of each node; the depth of color reflecting the rank in the hub genes.
About this article
Cite this article
Li, J., Shi, H., Yuan, Z. et al. The role of SPI1-TYROBP-FCER1G network in oncogenesis and prognosis of osteosarcoma, and its association with immune infiltration. BMC Cancer 22, 108 (2022). https://doi.org/10.1186/s12885-022-09216-w
- Immune Infiltration