Osteosarcoma (OS) is the most common primary malignant bone sarcoma, disproportionately affecting children and adolescents. It is a highly aggressive tumor: roughly 20% of patients already have detectable metastases at the time of initial diagnosis, and approximately 50% develop lung metastases at a later stage. Despite decades of multimodal treatment, the 5-year survival rate for patients with metastatic OS remains below 30%, a figure that has not substantially improved since the introduction of neoadjuvant chemotherapy in the 1980s. Identifying the molecular mechanisms that drive both the initial transformation and subsequent metastatic spread is therefore a clinical priority.
From single genes to networks: Prior research on OS molecular biology has largely focused on the role of individual oncogenes and tumor suppressors, including p53 (cell cycle arrest and apoptosis), SOX2 (stemness and migration), MALAT1 (proliferation and metastasis), and IGF-2 (osteoblast differentiation). While these studies have yielded important insights, single-gene analyses capture only a fraction of the complexity underlying OS behavior. More recent approaches have shifted toward identifying co-expressed gene clusters and protein-protein interaction (PPI) networks that collectively characterize OS oncogenesis and prognosis.
The immunological angle: A consistent observation across OS datasets is that the tumor immune microenvironment is substantially altered compared to normal bone tissue. Reduced immune cell infiltration and attenuated antigen presentation capacity are recurring features of OS, and both are linked to worse clinical outcomes. This study was designed to identify the specific gene network responsible for the immunological suppression observed in OS, and to determine whether that network also carries independent prognostic value.
Published in BMC Cancer in 2022, this study by Li, Shi, Yuan, and colleagues integrates multiple GEO microarray datasets with TCGA and TARGET survival data to identify a conserved immune-regulatory gene network centered on TYROBP, and demonstrates that its downregulation in OS correlates with attenuated immune infiltration and poorer overall survival.
The study drew on four GEO microarray datasets: GSE33382, GSE12865, GSE16088, and GSE14359. Each dataset was selected based on three inclusion criteria: it had to include both normal control and osteosarcoma tissue samples, come from Homo sapiens, and contain more than 10 samples. The control tissues differed meaningfully across datasets. In GSE12865 and GSE14359, controls were normal primary osteoblasts. In GSE33382, controls were osteoblasts derived from osteogenic differentiated bone-marrow-derived mesenchymal stem cells (MSCs). In GSE16088, the control samples came from an unspecified osteoblast cell line (H-012706). This heterogeneity was acknowledged as both a limitation and a strength, since overlapping results across heterogeneous controls are more likely to reflect true biological differences in OS rather than artifacts of any particular control cell type.
Differentially expressed gene (DEG) identification: GEO2R and the limma R package were used to identify DEGs in each dataset independently. Genes were considered significantly differentially expressed if they met a Benjamini-Hochberg (BH) false discovery rate (FDR) adjusted p-value below 0.05 and an absolute log2 fold change greater than 0.58. To reduce false positives from any individual dataset, only genes appearing as DEGs in three or more of the four datasets were carried forward into downstream analyses. This stringency yielded 124 upregulated genes and 98 downregulated genes.
PPI network and cluster analysis: STRING database was used to construct a PPI network from the full DEG set, which was then visualized and analyzed in Cytoscape. The MCODE plugin identified the top five densely interconnected clusters within the network. Hub genes were ranked by degree of connectivity, with the top 10 by degree calculated using the cytoHubba plugin with the Maximal Clique Centrality (MCC) algorithm. Survival analyses used the TARGET-OS dataset for primary analysis and GSE21257 for validation. Univariate Cox regression identified prognostically significant genes, and LASSO Cox regression (via the glmnet R package) was used to select an optimal multi-gene signature with penalized variable selection.
Immune infiltration estimation: Two complementary platforms were used to quantify immune cell composition in the tumor microenvironment. The xCell platform uses gene signatures to estimate proportions of 64 immune and stromal cell types, while ssGSEA (via the GSVA R package) calculates enrichment scores for immune cell gene sets at the sample level. GEPIA2 and the TCGA Pan-Cancer TPM dataset (from UCSC Xena Browser) were used to assess TYROBP co-expression across cancer types, and ChEA3 was used to predict transcription factors regulating the network.
After filtering for genes that were differentially expressed in at least three of the four GEO datasets, the authors identified 124 upregulated and 98 downregulated genes in OS compared to normal controls. Functional enrichment analysis of these two gene sets revealed a striking divergence in biological themes. The 124 upregulated genes were enriched in annotations related to collagen-containing extracellular matrix, smooth muscle cell proliferation, cell-substrate adhesion, and the p53 signaling pathway. These terms are consistent with the known mesenchymal origin of OS and its ability to remodel the surrounding stroma to facilitate invasion.
Downregulated genes and immune suppression: In contrast, the 98 downregulated genes were enriched in immune-related functions, including MHC class II protein complex and complex binding, lymphocyte proliferation, and leukocyte proliferation. This pattern indicated that OS samples systematically lose immunological surveillance capacity relative to normal osteoblast controls, with MHC class II antigen presentation capacity being particularly attenuated.
PPI network hub genes: When both upregulated and downregulated DEGs were combined into a PPI network in Cytoscape, the top 10 hub genes by connectivity degree were FN1, TYROBP, EGFR, CSF1R, C1QB, PLEK, LCP2, C1QA, CD86, and VWF. These hubs collectively showed enrichment in osteoclast differentiation, integrin binding, collagen-containing extracellular matrix, and glial cell activation. Notably, FN1 had the highest connectivity degree of all nodes but, when tested independently in survival analysis, did not significantly predict OS prognosis, suggesting that high connectivity does not automatically confer prognostic relevance.
Top 5 cluster functions: The five densely interconnected clusters had distinct biological identities. Top 1 cluster genes were all downregulated in OS and enriched in positive regulation of immune cell activation, cytokine production, and MHC class II protein complex. Top 2 cluster genes were associated with lipoprotein particle binding, granule membrane functions, and phagocytosis. Top 3 cluster was linked to actin cytoskeleton reorganization. Top 4 cluster was associated with fibroblast proliferation, and Top 5 cluster with smooth muscle cell migration and leukocyte migration.
To determine whether the functionally defined clusters also carried prognostic relevance, the authors performed univariate Cox regression across all 56 genes from the five top clusters using the TARGET-OS dataset. Nine of the 56 cluster genes were identified as significantly prognostic. Of these nine, five came from the TYROBP-centered Top 1 cluster: FCER1G, CSF1R, C1QA, TYROBP, and C1QB. All five had a hazard ratio (HR) below 1.0, meaning that higher expression of each gene was associated with better overall survival. This protective association aligns with the cluster's immune-activation function: higher expression means more immunological activity, which correlates with better outcomes.
TYROBP survival association: Kaplan-Meier analysis using median TYROBP expression as a cutoff showed that OS patients with higher TYROBP expression had significantly longer overall survival than those with lower expression. GSEA enrichment analysis further showed that genes positively co-expressed with TYROBP in OS were significantly enriched in innate immune system, adaptive immune system, cytokine signaling, neutrophil degranulation, and TYROBP causal network pathways from the Reactome gene set collection.
Cross-cancer TYROBP expression patterns: TYROBP expression was not uniformly altered across cancer types. Using TCGA Pan-Cancer TPM data from UCSC Xena, the authors found that TYROBP was significantly downregulated (adjusted p < 0.05) in OS, rectum adenocarcinoma (READ), pancreatic adenocarcinoma (PAAD), lung adenocarcinoma (LUAD), lung squamous cell carcinoma (LUSC), and colon adenocarcinoma (COAD). Conversely, TYROBP was upregulated in thyroid carcinoma (THCA), breast invasive carcinoma (BRCA), stomach adenocarcinoma (STAD), esophageal carcinoma (ESCA), and glioblastoma multiforme (GBM). This bidirectional pattern across cancer types suggested that the functional consequences of TYROBP expression depend on the specific tumor context and which immune cell types predominantly infiltrate each cancer type.
Conserved co-expression network: By intersecting the genes co-expressed with TYROBP (Spearman R > 0.6, p < 0.05) in cancer types where it is downregulated with those where it is upregulated, the authors identified 54 conserved TYROBP co-expression genes shared across all 11 examined cancer types. These 54 genes were enriched in mononuclear cell proliferation, lymphocyte proliferation, leukocyte proliferation, T cell activation, and antigen processing and presentation.
To build a parsimonious prognostic model from the 54-gene conserved TYROBP co-expression network, the authors first identified the top 10 hub genes within this network by PPI analysis using the cytoHubba MCC algorithm. The top 10 were SPI1 (first rank), TYROBP (second rank), FCER1G (third rank), ITGB2, C1QB, C1QA, LY86, LCP2, CCR1, and AIF1. The naming of the entire network as the "SPI1-TYROBP-FCER1G network" reflects the top three hub rankings, with TYROBP at the functional center as the signaling adaptor molecule.
Cross-cancer expression correlation: Across all 11 cancer types analyzed, the expression levels of SPI1, TYROBP, and FCER1G were highly correlated with each other, with pairwise Spearman R values exceeding 0.8 and p-values of 0 (below any computable threshold), confirming that these three genes reliably co-vary regardless of cancer type. This degree of co-expression across diverse tumor contexts strongly suggests they participate in a functionally coordinated regulatory module rather than being coincidentally co-expressed.
LASSO-derived risk score: LASSO Cox regression on the TARGET-OS dataset selected an optimal two-gene signature from the 10 hub genes. The final risk score formula was: (-0.1405) x FCER1G expression + (-0.0282) x SPI1 expression. Both coefficients are negative, meaning that higher expression of either gene reduces the risk score, and lower risk scores correspond to better survival outcomes. The AUC of this combined risk score for overall survival was 0.673, which was higher than the individual AUCs of FCER1G alone (0.664), TYROBP alone (0.634), and SPI1 alone (0.658), demonstrating additive prognostic value from combining the two genes.
Kaplan-Meier validation: Dividing patients by median risk score showed that those with higher risk scores (lower SPI1 and FCER1G expression) had significantly shorter overall survival. Although the AUC of 0.673 is modest in absolute terms, it is derived from a biologically coherent two-gene model applied to a pediatric sarcoma dataset with inherently limited sample sizes, and it outperforms either gene alone.
Beyond the continuous risk score, the authors used consensus clustering on SPI1 and FCER1G expression to divide OS patients into two discrete clusters, a strategy that avoids the arbitrariness of median cutoffs and allows discovery of natural groupings in the data. In the primary TARGET-OS dataset, Cluster 1 (lower SPI1 and FCER1G expression) was associated with significantly poorer overall survival (p = 0.010) and significantly lower stromal scores (p < 0.01) and immune scores (p < 0.001) than Cluster 2.
Immune cell composition in TARGET-OS: Using the ssGSEA algorithm, almost all immune cell types showed higher enrichment scores in Cluster 2 than Cluster 1, with the exceptions of CD56dim NK cells and Th2 cells. The xCell platform similarly found higher proportions of CD4+ effector memory T cells (Tem), endothelial cells (both lymphatic and microvascular), M1 and M2 macrophages, monocytes, and multiple dendritic cell subsets (activated DC, conventional DC, immature DC, plasmacytoid DC) in Cluster 2, after excluding cell types with proportions below 1%.
Validation in GSE21257: The same clustering approach was applied to the independent GSE21257 dataset, and the results were consistent. Cluster 1 again had significantly poorer overall survival (p = 0.021) and lower ssGSEA immune scores for nearly all cell types. xCell identified higher immune scores and microenvironment scores but lower stroma scores in Cluster 2, confirming that the SPI1/FCER1G-low phenotype represents a less immunologically active, more stromal tumor microenvironment.
Metastasis correlation: Supporting the clinical relevance of this network, comparison of OS patients with and without metastases in GSE21257 showed that all three hub genes were significantly downregulated in metastatic tumors: TYROBP (p = 0.0003, Log2FC = -1.16), FCER1G (p = 0.0002, Log2FC = -1.19), and SPI1 (p = 0.0059, Log2FC = -0.599). The stronger downregulation of TYROBP and FCER1G compared to SPI1 in metastatic disease is consistent with the known role of ITAM-containing adaptor proteins in mediating immune cell recruitment at distant sites.
Although TYROBP sits at the center of the co-expression network by PPI hub ranking, it cannot directly regulate gene expression because it lacks a DNA-binding domain. To identify the transcription factor initiating the network, the authors used the ChEA3 platform, which performs transcription factor enrichment analysis using orthogonal omics data from multiple databases. SPI1 (also known as PU.1) ranked second among the top 50 predicted transcription factors for the 55-gene TYROBP co-expression network. Critically, SPI1 was found to bind all but five genes in the network, including TYROBP and FCER1G themselves, making SPI1 a plausible master regulator of the entire module.
ChIP-seq validation: ENCODE ChIP-seq data from multiple cell lines (GM12878, GM12891, K562, HL-60) and ChIP-seq data from THP-1 cells (a monocytic cell line) confirmed direct binding of SPI1/PU.1 to the TYROBP promoter region and at least nine other network hub gene loci. GSEA enrichment analysis independently showed that PU.1 and IRF transcription factor binding motifs were significantly enriched in genes positively co-expressed with TYROBP in OS. Low SPI1 expression in the TARGET-OS dataset was significantly associated with poorer overall survival, further supporting its functional role.
Upstream regulation by TNF-alpha: GSEA on TYROBP co-expressed genes showed significant enrichment of "TNF-alpha signaling via NF-kB," "IFN-alpha/gamma response," "inflammatory response," and "IL2-STAT5 signaling" hallmark gene sets, suggesting that the network is activated by inflammatory cytokine signaling in the tumor microenvironment. Drug Pair Seeker (DPS) analysis identified TNF-alpha and IGF2 as cytokine pairs capable of partially rescuing the underexpression of the TYROBP network in OS. This finding has an important clinical implication: anti-TNF-alpha therapies, which are used for inflammatory conditions, may worsen OS prognosis by further suppressing SPI1-TYROBP-FCER1G network activity.
Downstream immune effectors: Overlapping TYROBP positively co-expressed genes (R > 0.4, p < 0.05) with OS-downregulated genes across 2 or more GEO datasets identified specific immune effector molecules suppressed by network downregulation. These included one immune stimulator (CD86), five MHC molecules (HLA-DMA, HLA-DMB, HLA-DPA1, HLA-DPB1, HLA-DRA), three chemokines (CCL4, CXCL10, CX3CL1), and two immune inhibitors (CSF1R and HAVCR2/TIM-3). Loss of chemokines CCL4 (recruiting CD8+ T cells and NK cells via CCR5) and CXCL10 (recruiting CD8+ T cells and NK cells via CXCR3) reduces immune cell trafficking to the tumor site. Loss of MHC class II molecules impairs T cell activation once immune cells arrive.
Dataset heterogeneity and sample size constraints: The four GEO microarray datasets used in DEG analysis differed in the nature of their control samples, ranging from primary osteoblasts to MSC-derived osteoblasts to a cell line, and the OS samples themselves included both primary biopsies and metastatic tissues. The authors acknowledge that this heterogeneity likely reduced the number of overlapping DEGs identified across datasets. Using a three-dataset overlap threshold was a deliberate conservative choice that prioritizes reproducibility at the cost of sensitivity, meaning some genuine OS-associated genes may have been excluded. The TARGET-OS and GSE21257 survival datasets, while the largest publicly available OS survival cohorts, are still relatively small by oncology standards, which limits the statistical power of survival modeling and increases the risk that the AUC of 0.673 may not generalize to larger prospective cohorts.
Bioinformatic findings require experimental validation: The entire study is computational. The proposed mechanistic pathway, in which TNF-alpha activates SPI1, which in turn transcriptionally drives TYROBP, FCER1G, and the downstream immune effectors, is strongly supported by bioinformatics but has not yet been confirmed by in vitro or in vivo experiments. Key unknowns include whether the SPI1-TYROBP-FCER1G network is expressed in OS tumor cells themselves, in the myeloid cells and endothelial cells of the tumor microenvironment, or in both. The relative contributions of each cell compartment to the network's prognostic signal remain uncharacterized.
Cancer-type-specific directionality: One of the most intriguing and unresolved observations in this paper is that the SPI1-TYROBP-FCER1G network is downregulated in OS but upregulated in cancers such as breast cancer, clear cell renal cell carcinoma, low-grade glioma, and gastric cancer. In those cancers, higher network expression is associated with higher immune infiltration but, paradoxically, with worse prognosis in some contexts. The authors cannot explain from their data why the same network has opposing prognostic implications across cancer types, and suggest this may reflect differences in the dominant infiltrating immune cell types or differences in how immune infiltration affects tumor progression in each context.
Therapeutic targeting potential: The DPS analysis suggesting TNF-alpha and IGF2 as potential rescuers of the suppressed TYROBP network provides a starting point for therapeutic hypotheses. TYROBP and FCER1G both contain immunoreceptor tyrosine-based activation motifs (ITAM), making them potentially targetable by agents that modulate ITAM-mediated signaling in myeloid and NK cells. The downregulation of CSF1R and TIM-3 (HAVCR2) in Cluster 1 tumors also points toward the tumor microenvironment as a therapeutic target, though their suppression in OS differs mechanistically from the immune checkpoint context in other cancers. Future work combining single-cell RNA sequencing with functional perturbation experiments will be needed to confirm the causal role of this network in OS immune escape and metastasis.