Skip to main content

The role of SPI1-TYROBP-FCER1G network in oncogenesis and prognosis of osteosarcoma, and its association with immune infiltration

Abstract

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.

Peer Review reports

Introduction

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 [3]. Now the treatment for OS has reached a plateau to the patients with metastases, the 5-year survival rate of which is under 30% [4]. 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) [5], SOX2 (enhancing cell stemness and migration) [6], MALAT1 (promoting proliferation and metastasis) [7], IGF-2 (affecting osteoblast differentiation) [8], cyclin E1 (inhibiting proliferation and increasing chemotherapeutics sensitivity) [9], and PTH (affecting migration) [10]. 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 [16]) 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 [17].

ENCODE PU.1 ChIP-seq peaks data from multiple cell lines (i.e., GM12878, GM12891, K562, HL-60) were downloaded [18] and visualized by software IGV [19]. 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 [20]. 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 [21].

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 [24].

Functional enrichment analysis

Gene Ontology Analysis on Gene Clusters Gene ontology (GO) (including biological process (BP), molecular functions (MF), and cellular components (CC) terms) [25] and Kyoto Encyclopedia of Genes and Genomes (KEGG) pathways [26] are used to annotate the functions of genes via “clusterProfiler” R package [27] 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 [28] 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 [29]. The network data of proteins from DEGs were downloaded from STRING to visualize the protein interaction network in the Cytoscape software [30]. 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 [31]. 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) [32] in Cytoscape.

Survival analyses

Survival analyses were performed using univariate Cox regression and LASSO regression via “survival” [33] and “glmnet” [34] 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” [35] and “ggplot2” R package.

Immune infiltration analysis

Platform xCell (xCell signature) [36] and single-sample gene set enrichment analysis (ssGSEA, based on immune signature by “gsva” R package [37]) 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.

Results

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).

Fig. 1
figure 1

Differentially expressed genes (DEGs) were identified in four OS datasets and functional enrichment analysis on them. A-B The Venn diagram representing DEGs numbers in each dataset. C-D The top three KEGG pathways and GO functions (BP, CC, MF) assembled by upregulated (C) and downregulated (D) genes

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.

Fig. 2
figure 2

The top 10 hub genes in PPI network of DEGs and the top 5 cluster genes in DEGs network. A The top 10 hub genes in 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. B The top 3 KEGG and GO annotations enriched by the total DEGs. C-G The PPI network of top 5 clusters

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).

Fig. 3
figure 3

The prognostic genes in the top 5 clusters. A The Venn diagram representing the 9 overlapping genes between 56 of top 5 cluster genes and all prognostic genes identified by univariate Cox regression analysis on TARGET-OS dataset. B The forest plot showing the Hazard Ratio with Credible Interval and P value of each gene in the analysis. C Survival curve showing OS overall survival based on the median expressional level cutoff of TYROBP. D GSEA plot showing significant enrichment of immune-related Reactome pathways in TYROBP positively co-expressed genes. E Survival curve showing OS overall survival based on the median expressional level cutoff of FN1

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.

Fig. 4
figure 4

The common 55 TYROBP co-expression genes (including TYROBP) across cancer types and their enriched functions. A The boxplot showing the expressional levels of TYROBP between para-cancer tissues and cancer tissues across different cancer types (*, p < 0.05; **, p < 0.01; ***, p < 0.001). B The Venn diagram representing the number of TYROBP co-expression genes in cancer types with upregulated TYROBP (Left) and cancer type with downregulated TYROBP (Right), and the 55 of common genes across all cancer types. C The top 20 in all annotations and the top 3 in each type of database (BP, CC, MF, KEGG) enriched by the 55 TYROBP co-expression genes

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).

Fig. 5
figure 5

The prognostic two-gene signature from the hub genes in conserved SPI1-TYROBP-FCER1G network. A The hub genes in conserved SPI1-TYROBP-FCER1G network. The depth of color reflecting the rank in the hub genes. B The correlation between SPI1/FCER1G and TYROBP across cancer types. C Coefficients of selected features are shown by lambda parameter. Partial likelihood deviance versus log (λ) was performed through LASSO regression. D ROC curve based on the two-gene signature for overall survival according to risk scores. E ROC curves for overall survival according to the expressional levels of SPI1, FCER1G, or TYROBP. F Kaplan–Meier plots of overall survival according to risk scores

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).

Fig. 6
figure 6

Differential overall survival and immune infiltration of OS in cluster 1 and cluster 2 subgroups according to the expressional levels of SPI1 and FCER1G. A-E Based on TARGET-OS dataset. F-J Based on GSE21257 dataset. A and F Consensus clustering matrix for k = 2. B and G Kaplan–Meier overall survival curves for OS. C and H The boxplot of immune cell infiltration score in cluster 1 and 2. D and I The boxplot of immune cell infiltration proportion in cluster 1 and 2. E and J The boxplot of immune infiltration score in cluster 1 and 2. ns, p ≥ 0.05,*, p < 0.05,**, p < 0.01,***, p < 0.001

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).

Fig. 7
figure 7

The key transcription factor SPI1 in the conserved SPI1-TYROBP-FCER1G network. A The Venn diagram representing the numbers of ChEA3 predicted transcription factors/SPI1 binding genes and TYROBP network genes. B GSEA plot showing significant enrichment of PU.1 and IRF binding motifs in TYROBP positively co-expressed genes. C The genomic signal density around TYROBP showing the PU.1 binding regions upstream and the mRNA level of TYROBP. D Kaplan–Meier plots of overall survival according to expressional level of SPI1 (median cutoff). E GSEA plot showing the significant enrichments of signaling pathways including TNF-α signaling via NF-κB in TYROBP positively co-expressed genes. F Drug pairs found by Drug Pair Seeker representing the cytokine pairs including TNF-α that can upregulate the TYROPB network

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).

Fig. 8
figure 8

The relationship between TYROBP network and immune-related genes in OS and the schematic workflow of the whole study. A-B The Venn diagram representing the overlapping between downregulated genes in OS, TYROBP positively correlated genes in OS, and immune related genes (i.e., immune stimulators, MHC molecules, chemokines in A, and immune inhibitors in B). C The working model representing the TNF-α-SPI1-TYROBP network pathway across cancer types, which can regulate immune infiltration and CD8+ T cell activation. D Workflow for the whole study as to how to identify the SPI1-TYROBP-FCER1G network and its association with immune microenvironment

Discussion

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 [52]), molecules CD86 (co-stimulatory molecule CD86 polymorphism is also associated with increased susceptibility to osteosarcoma [53]) 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 [68] 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).

Abbreviations

aDC:

Activated dendritic cells

cDC:

Conventional dendritic cells

iDC:

Immature dendritic cells

pDC:

Plasmacytoid dendritic cells

Tcm:

Central memory T-cells

Tem:

Effector memory T cells

Tgd:

T cells gamma delta

CMP:

Common myeloid progenitors

GMP:

Granulocyte–macrophage progenitors

MSC:

Mesenchymal stem cells

CLP:

Common lymphoid progenitors

ly endothelial:

Lymphatic endothelial cells

mv endothelial:

Microvascular endothelial cells

OS:

Osteosarcoma

THCA:

Thyroid carcinoma

BRCA:

Breast invasive carcinoma

ESCA:

Esophageal carcinoma

GBM:

Glioblastoma multiforme

STAD:

Stomach adenocarcinoma

PAAD:

Pancreatic adenocarcinoma

READ:

Rectum adenocarcinoma

COAD:

Colon adenocarcinoma

LUAD:

Lung adenocarcinoma

LUSC:

Lung squamous cell carcinoma

References

  1. Raymond AK, Jaffe N. Osteosarcoma multidisciplinary approach to the management from the pathologist’s perspective. Cancer Treat Res. 2009;152:63–84.

  2. Zhao X, et al. Osteosarcoma: a review of current and future therapeutic approaches. Biomed Eng Online. 2021;20(1):24.

    Article  PubMed  PubMed Central  Google Scholar 

  3. Tang N, et al. Osteosarcoma Development and Stem Cell Differentiation. Clin Orthop Relat Res. 2008;466:2114–30.

    Article  PubMed  PubMed Central  Google Scholar 

  4. Ferguson JL, Turner SP. Bone Cancer: Diagnosis and Treatment Principles. Am Fam Physician. 2018;98(4):205–13.

    PubMed  Google Scholar 

  5. Chen Z, et al. TP53 Mutations and Survival in Osteosarcoma Patients: A Meta-Analysis of Published Data. Dis Markers. 2016;2016:4639575.

    PubMed  PubMed Central  Google Scholar 

  6. Maurizi G, et al. Sox2 is required for tumor development and cancer cell proliferation in osteosarcoma. Oncogene. 2018;37:4626–32.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  7. 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.

    CAS  PubMed  Google Scholar 

  8. 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.

    Article  CAS  PubMed  Google Scholar 

  9. Wei R, et al. Cyclin E1 is a prognostic biomarker and potential therapeutic target in osteosarcoma. J Orthop Res. 2020;38(9):1952–64.

    Article  CAS  PubMed  Google Scholar 

  10. 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.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  11. Jin Z, et al. Cross-Species Gene Expression Analysis Reveals Gene Modules Implicated in Human Osteosarcoma. Front Genet. 2019;10:697.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  12. Niu J, et al. Identification of Potential Therapeutic Targets and Immune Cell Infiltration Characteristics in Osteosarcoma Using Bioinformatics Strategy. Front Oncol. 2020;10:1628.

    Article  PubMed  PubMed Central  Google Scholar 

  13. 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.

    Article  CAS  PubMed  Google Scholar 

  14. Jia Y, et al. Identification of potential gene signatures associated with osteosarcoma by integrated bioinformatics analysis. PeerJ. 2021;9:e11496.

    Article  PubMed  PubMed Central  Google Scholar 

  15. Liu J, et al. Identification of potential crucial genes and key pathways in osteosarcoma. Hereditas. 2020;157(1):29.

    Article  PubMed  PubMed Central  Google Scholar 

  16. Vivian J, et al. Toil enables reproducible, open source, big biomedical data analyses. Nat Biotechnol. 2017;35(4):314–6.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  17. 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.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  18. Davis CA, et al. The Encyclopedia of DNA elements (ENCODE): data portal update. Nucleic Acids Res. 2018;46(D1):D794-d801.

    Article  CAS  PubMed  Google Scholar 

  19. Robinson JT, et al. Integrative genomics viewer. Nat Biotechnol. 2011;29(1):24–6.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  20. 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.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  21. Keenan AB, et al. ChEA3: transcription factor enrichment analysis by orthogonal omics integration. Nucleic Acids Res. 2019;47(W1):W212–24.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  22. Davis S, Meltzer PS. GEOquery: a bridge between the Gene Expression Omnibus (GEO) and BioConductor. Bioinformatics. 2007;23(14):1846–7.

    Article  PubMed  CAS  Google Scholar 

  23. Ritchie ME, et al. limma powers differential expression analyses for RNA-sequencing and microarray studies. Nucleic Acids Res. 2015;43(7):e47–e47.

    Article  PubMed  PubMed Central  CAS  Google Scholar 

  24. Gómez-Rubio V. ggplot2-elegant graphics for data analysis. J Stat Softw. 2017;77:1–3.

    Article  Google Scholar 

  25. Ashburner M, et al. Gene ontology: tool for the unification of biology The Gene Ontology Consortium. Nat Genet. 2000;25(1):25–9.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  26. Kanehisa M, Goto S. KEGG: kyoto encyclopedia of genes and genomes. Nucleic Acids Res. 2000;28(1):27–30.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  27. Yu G, et al. clusterProfiler: an R package for comparing biological themes among gene clusters. Omics. 2012;16(5):284–7.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  28. 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.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  29. 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.

    Article  CAS  PubMed  Google Scholar 

  30. Shannon P, et al. Cytoscape: a software environment for integrated models of biomolecular interaction networks. Genome Res. 2003;13(11):2498–504.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  31. Bader GD, Hogue CWV. An automated method for finding molecular complexes in large protein interaction networks. BMC Bioinformatics. 2003;4(1):2.

    Article  PubMed  PubMed Central  Google Scholar 

  32. 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.

    Article  PubMed  Google Scholar 

  33. Therneau TM, Lumley T. Package ‘survival.’ R Top Doc. 2015;128(10):28–33.

  34. Simon N, et al. Regularization Paths for Cox’s Proportional Hazards Model via Coordinate Descent. J Stat Softw. 2011;39(5):1–13.

  35. Robin X, et al. pROC: An open-source package for R and S+ to analyze and compare ROC curves. BMC Bioinformatics. 2011;12:77.

    Article  PubMed  PubMed Central  Google Scholar 

  36. Aran D, Hu Z, Butte AJ. xCell: digitally portraying the tissue cellular heterogeneity landscape. Genome Biol. 2017;18(1):220.

    Article  PubMed  PubMed Central  CAS  Google Scholar 

  37. Hänzelmann S, Castelo R, Guinney J. GSVA: gene set variation analysis for microarray and RNA-seq data. BMC Bioinformatics. 2013;14:7.

    Article  PubMed  PubMed Central  Google Scholar 

  38. 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.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  39. 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.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  40. 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.

    Article  CAS  PubMed  Google Scholar 

  41. Wu P, et al. TYROBP is a potential prognostic biomarker of clear cell renal cell carcinoma. FEBS Open Bio. 2020;10(12):2588–604.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  42. 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.

    Article  PubMed  PubMed Central  Google Scholar 

  43. 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.

    Article  PubMed  PubMed Central  CAS  Google Scholar 

  44. Liang T, et al. TYROBP, TLR4 and ITGAM regulated macrophages polarization and immune checkpoints expression in osteosarcoma. Sci Rep. 2021;11(1):19315.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  45. Chen H, et al. Three-gene prognostic biomarkers for seminoma identified by weighted gene co-expression network analysis. PLoS One. 2020;15(10):e0240943.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  46. Mukaida N, Sasaki SI, Baba T. CCL4 Signaling in the Tumor Microenvironment. Adv Exp Med Biol. 2020;1231:23–32.

    Article  CAS  PubMed  Google Scholar 

  47. Yang H, et al. A comprehensive analysis of immune infiltration in the tumor microenvironment of osteosarcoma. Cancer Med. 2021;10(16):5696–711.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  48. Flores RJ, et al. A novel prognostic model for osteosarcoma using circulating CXCL10 and FLT3LG. Cancer. 2017;123(1):144–54.

    Article  CAS  PubMed  Google Scholar 

  49. Pradelli E, et al. Antagonism of chemokine receptor CXCR3 inhibits osteosarcoma metastasis to lungs. Int J Cancer. 2009;125(11):2586–94.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  50. 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.

    Article  CAS  PubMed  Google Scholar 

  51. Conroy MJ, Lysaght J. CX3CL1 signaling in the tumor microenvironment. Tumor Microenvironment. 2020;1231:1–12.

    Article  CAS  Google Scholar 

  52. 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.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  53. Wang W, et al. CD86 + 1057G/A polymorphism and susceptibility to osteosarcoma. DNA Cell Biol. 2011;30(11):925–9.

    Article  CAS  PubMed  Google Scholar 

  54. 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.

    CAS  PubMed  PubMed Central  Google Scholar 

  55. 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.

    Article  CAS  PubMed  Google Scholar 

  56. Anderson AC. Tim-3, a negative regulator of anti-tumor immunity. Curr Opin Immunol. 2012;24(2):213–6.

    Article  CAS  PubMed  Google Scholar 

  57. Ge W, et al. Tim-3 as a diagnostic and prognostic biomarker of osteosarcoma. Tumour Biol. 2017;39(7):1010428317715643.

    Article  PubMed  CAS  Google Scholar 

  58. 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.

    Article  CAS  PubMed  Google Scholar 

  59. 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.

    Article  CAS  Google Scholar 

  60. Dang D, et al. Computational Approach to Identifying Universal Macrophage Biomarkers. Front Physiol. 2020;11:275.

    Article  PubMed  PubMed Central  Google Scholar 

  61. Pan Y-G, et al. FcεRI γ-Chain Negatively Modulates Dectin-1 Responses in Dendritic Cells. Frontiers in Immunology. 2017;8:1424.

    Article  PubMed  PubMed Central  CAS  Google Scholar 

  62. Hamerman JA, et al. The expanding roles of ITAM adapters FcRgamma and DAP12 in myeloid cells. Immunol Rev. 2009;232(1):42–58.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  63. 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.

    Article  CAS  PubMed  Google Scholar 

  64. Nutt SL. Directing the conductor: TNF regulation of HSCs. Blood, The Journal of the American Society of Hematology. 2019;133(8):771–3.

    CAS  Google Scholar 

  65. 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.

    Article  CAS  PubMed  Google Scholar 

  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.

    CAS  PubMed  Google Scholar 

  67. 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.

    Article  CAS  PubMed  Google Scholar 

  68. Kirkham B, et al. Tumor necrosis factor-alpha inhibitors: An overview of adverse effects. Waltham, MA: UpToDate; 2016.

    Google Scholar 

Download references

Acknowledgements

Not applicable.

Funding

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).

Author information

Authors and Affiliations

Authors

Contributions

JL, HS, ZY and ML2: searched data, analyzed data, and wrote the manuscript. ZW, HL, YL: analyzed data. ML3 and ML2: edited the manuscript. ML3 and ML2: conception of the project. The author(s) read and approved the final manuscript.

Corresponding authors

Correspondence to Ming Lu or Ming Lu.

Ethics declarations

Ethics approval and consent to participate

All methods were carried out in accordance with relevant guidelines and regulations (Declaration of Helsinki).

Consent for publication

Not applicable.

Competing interests

The authors declare that they have no competing interests.

Additional information

Publisher's Note

Springer Nature remains neutral with regard to jurisdictional claims in published maps and institutional affiliations.

Supplementary Information

Additional file 1: Figure 1.

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.

Rights and permissions

Open Access This article is licensed under a Creative Commons Attribution 4.0 International License, which permits use, sharing, adaptation, distribution and reproduction in any medium or format, as long as you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons licence, and indicate if changes were made. The images or other third party material in this article are included in the article's Creative Commons licence, unless indicated otherwise in a credit line to the material. If material is not included in the article's Creative Commons licence and your intended use is not permitted by statutory regulation or exceeds the permitted use, you will need to obtain permission directly from the copyright holder. To view a copy of this licence, visit http://creativecommons.org/licenses/by/4.0/. The Creative Commons Public Domain Dedication waiver (http://creativecommons.org/publicdomain/zero/1.0/) applies to the data made available in this article, unless otherwise stated in a credit line to the data.

Reprints and permissions

About this article

Check for updates. Verify currency and authenticity via CrossMark

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

Download citation

  • Received:

  • Accepted:

  • Published:

  • DOI: https://doi.org/10.1186/s12885-022-09216-w

Keywords