Osteosarcoma is the most prevalent primary bone malignancy and the 8th most frequent cancer overall, disproportionately striking children and young adults. While combining surgery with chemotherapy has pushed long-term survival to roughly 60-70% in patients without distant metastasis, that figure drops precipitously once metastatic spread occurs, and the underlying causes of why some tumors stay confined while others disseminate remain poorly understood. This knowledge gap prevents effective early interception and discourages the development of targeted therapies.
The core challenge is biological heterogeneity. Osteosarcoma arises from mesenchymal stem cells and accumulates a chaotic landscape of genomic rearrangements rather than a tidy set of recurrent driver mutations. This makes it hard to define a single molecular signature that reliably separates normal bone from primary tumor, or primary tumor from metastatic lesion. Identifying genes that are consistently deregulated across multiple independent patient cohorts would give clinicians and researchers stable targets for biomarker development and drug design.
Bioinformatics offers a path forward. The Gene Expression Omnibus (GEO) database aggregates microarray and RNA-seq datasets from laboratories around the world, and integrating multiple osteosarcoma datasets in a single analysis allows researchers to filter out study-specific noise and focus on the signal that reproduces across platforms. This paper by Liu et al. (2020) applies that strategy, combining three GEO datasets totaling 112 osteosarcoma tissue samples and 7 osteoblast controls to systematically map differentially expressed genes (DEGs) at two critical transitions: normal-to-primary tumor and primary-to-metastatic tumor.
The authors applied a four-criterion filter to select appropriate GEO datasets. Eligible datasets had to include primary or metastatic osteosarcoma tissues as tumor samples, normal human bone or osteoblasts as comparators, at least 1,000 DEGs meeting FDR < 0.05 and |log2 fold-change| > 1, and at least 10 overlapping DEGs with the other chosen datasets. Three datasets passed all criteria: GSE14359 (Affymetrix HG U133A platform, Germany, comprising 2 osteoblasts, 10 primary tissue samples, and 4 lung metastases), GSE16088 (Affymetrix HG U133A, USA, with 3 osteoblast cell lines and 14 primary tissue samples), and GSE33383 (Illumina human-6 v2.0, Norway, with 3 osteoblasts and 84 primary tissue samples). The combined sample pool of 112 tumor tissues and 7 osteoblast controls is a meaningful size for a bioinformatics meta-analysis in a rare cancer.
Data preprocessing: Raw probe-level CEL files were processed using the Robust Multiarray Average (RMA) algorithm implemented in the R Affy package. RMA handles background correction and quantile normalization to place all samples on a common expression scale before any comparative statistics are run. Where a gene mapped to multiple probe sets, the values were averaged to produce a single representative expression value per gene per sample.
DEG identification: Differential expression was computed using the LIMMA (Linear Models for Microarray Data) package in R, which applies empirical Bayes moderation to improve statistical power with smaller sample sizes. The Benjamini-Hochberg method was used to control the false discovery rate (FDR), and the study set cut-offs at FDR < 0.05 combined with |log2FC| > 1, meaning only genes at least 2-fold different in expression and statistically robust were retained. Volcano plots visualized the full DEG landscape for each dataset, and Venn diagrams identified the overlapping genes across datasets, which are the reproducible candidates of greatest interest.
Downstream analyses: Reproducible DEGs were subjected to Gene Ontology (GO) enrichment using the clusterProfiler R package, covering biological process, cellular component, and molecular function categories, alongside KEGG pathway analysis. All GO and KEGG terms with FDR < 0.05 were considered significant. Protein-protein interaction (PPI) networks were built using the STRING database at a confidence score cutoff of 0.7, visualized in Cytoscape 3.7.2, and modular hubs were extracted using the MCODE (Molecular Complex Detection) plug-in to identify the most densely connected gene clusters.
Comparing osteoblast controls to primary osteosarcoma samples in each dataset individually produced large lists of normal-primary DEGs (NPDEGs): 777 in GSE14359 (489 up-regulated, 288 down-regulated), 1,943 in GSE16088 (1,010 up, 933 down), and 771 in GSE33383 (350 up, 421 down). The large variation in total DEG count across datasets reflects differences in platform technology, sample size, and the particular osteoblast cell lines used as normal references, which is exactly why cross-dataset filtering is essential. After applying Venn diagram overlap analysis, only 74 NPDEGs appeared in all three datasets simultaneously. These 74 genes represent the most consistent molecular signature of osteosarcoma tumorigenesis detectable across independent international cohorts.
Biological process enrichment: GO analysis of the 74 NPDEGs identified 82 significant biological process categories. The top enriched categories clustered around antigen processing and presentation, phagocytosis, and immune cell activation, particularly involving MHC class II receptors. This is biologically meaningful because the bone microenvironment is a specialized immune compartment where osteoclasts, which function as highly specialized macrophages, are tightly regulated by immunological cytokines. Disruption of immune surveillance in this niche appears to be a central feature of early osteosarcoma oncogenesis.
KEGG pathway enrichment: KEGG analysis reinforced the immune theme, with significant enrichment in pathways related to antigen processing and presentation via MHC class II molecules. This suggests that loss of normal immune recognition in the bone microenvironment, potentially impeding cytotoxic T-cell responses against nascent tumor cells, contributes to early osteosarcoma establishment. The 74 NPDEGs also included genes involved in extracellular matrix organization and cell-surface receptor signaling, reflecting the structural remodeling that accompanies tumor formation in bone.
The metastasis analysis was confined to GSE14359, the only dataset containing both primary and metastatic (lung) osteosarcoma samples alongside the normal osteoblast controls. Comparing primary tumors to metastatic lesions produced 764 primary-metastatic DEGs (PMDEGs): 309 up-regulated and 455 down-regulated. The larger number of down-regulated genes in this comparison is notable, potentially reflecting loss of tissue-specific programs in osteosarcoma cells that have successfully colonized the lung.
Biological process enrichment: GO analysis of the 764 PMDEGs identified 162 significant biological process categories. The single most significant category was mitotic nuclear division, reflecting increased proliferative activity in metastatic lesions. The top cellular component category was the extracellular matrix (ECM), which is consistent with the well-established requirement for ECM remodeling during tumor cell invasion and dissemination. The molecular function categories pointed toward cytoskeletal binding and motor activity, further supporting a migratory, invasive phenotype.
KEGG pathway enrichment: Three KEGG pathways stood out as most enriched in the PMDEGs: cell adhesion molecules (CAMs), focal adhesion, and ECM-receptor interaction. All three are canonical components of the metastatic cascade, governing how tumor cells detach from the primary site, survive in circulation, and re-anchor at distant organs. Additionally, the tumor necrosis factor (TNF) signaling pathway was significantly correlated with osteosarcoma metastasis. The TNF pathway is constitutively activated in human osteosarcoma cells and has been proposed as a chemotherapy target for advanced disease, making its appearance here clinically relevant.
The most clinically informative subset of findings came from comparing the up-regulated NPDEGs with the up-regulated PMDEGs using a Venn diagram. Seven genes appeared in both sets, meaning they were progressively up-regulated across all three disease stages: normal osteoblast, primary tumor, and metastatic lesion. These genes are VAMP8, A2M, HLA-DRA, SPARCL1, HLA-DQA1, APOC1, and AQP1. No genes were identified as continuously down-regulated across both transitions, which underscores the specificity of this ascending pattern. Genes that behave this way are particularly attractive as biomarkers because a single assay could, in principle, stratify patients along the full malignancy spectrum.
HLA-DRA: Of the seven, HLA-DRA (Major Histocompatibility Complex Class II, DR alpha) emerged as the most prominent candidate, being the hub node in the NPDEG PPI network module with the highest degree of connectivity. HLA-DRA encodes an MHC class II antigen-presenting molecule and is involved in 18 distinct GO biological process categories and multiple KEGG pathways including autoimmune disease, allograft rejection, and immune system-related signaling. Its sustained up-regulation during osteosarcoma progression may reflect an attempt by the tumor microenvironment to engage immune cells, or alternatively a tumor-cell-intrinsic hijacking of immune presentation machinery. Prior studies have established HLA-DRA as a predictor of osteosarcoma metastasis, lending external validation to its identification here.
SPARCL1: SPARCL1 (SPARC-Like 1) is an ECM remodeling gene that modulates extracellular calcium by binding to collagen I, pointing to a potential role in osteosarcoma cell metastasis through matrix remodeling. An interesting discrepancy exists in the literature: Zhao et al. reported that SPARCL1 is down-regulated in osteosarcoma via epigenetic methylation and that it suppresses metastasis, while this dataset shows SPARCL1 to be continuously up-regulated. The authors acknowledge this contradiction without resolution, suggesting that the expression behavior of SPARCL1 may be context-dependent or cohort-specific, and that further experimental validation is needed.
AQP1, APOC1, VAMP8, HLA-DQA1, A2M: The remaining five genes have diverse functions spanning vesicle trafficking (VAMP8), lipoprotein metabolism (APOC1), water channel biology (AQP1), immune antigen presentation (HLA-DQA1), and protease inhibition (A2M). Their shared upward trajectory across tumor stages suggests the osteosarcoma transcriptome activates a coordinated program involving immune evasion, lipid metabolism, and membrane trafficking as tumors progress.
Protein-protein interaction networks translate gene lists into functional maps by connecting proteins known to physically interact or functionally associate. The STRING database at a confidence cutoff of 0.7 was used to construct two separate networks: one from the 74 NPDEGs and one from the 764 PMDEGs. The NPDEG network contained 49 nodes and 91 interactions, a relatively sparse but coherent network given the small input gene set. MCODE analysis identified one significant module with a score of 5, and HLA-DRA was the hub with the highest degree of connectivity. GO enrichment of this module's genes with the ClueGO tool confirmed that the dominant biological process was antigen processing and presentation of exogenous peptide antigens via MHC class II, validating the immune biology interpretation from the pathway-level analysis.
PMDEG network hubs: The PMDEG network was far more complex, with 521 nodes and 2,415 interactions among the 764 PMDEGs. MCODE identified three significant modules, all with scores at or above 10. Module 1 (score = 32.5) included 36 genes and had CDK1, CDK20, and CCNB1 as hub nodes; its dominant GO category was nuclear division. Module 2 (score = 13.8) included 14 genes with MTIF2 and MRPS7 as hubs; its GO enrichment pointed to mitochondrial translation and energy metabolism. Module 3 included VEGFA and EGF as key nodes, both of which are well-known growth factor signaling drivers with established roles in tumor angiogenesis and proliferation.
Cell cycle hubs CDK1, CDK20, CCNB1: Cyclin-dependent kinase 1 (CDK1) and CDK20 belong to the serine/threonine protein kinase family, and Cyclin B1 (CCNB1) is a pivotal regulator of the G2/M cell cycle transition. All three are involved in cell cycle control and growth. Reduction of CDK1 activity has been shown to be critical for osteosarcoma cell survival, making CDK1 and its regulatory partners candidate therapeutic targets. VEGFA and EGF are the hub genes in the angiogenesis-related module, consistent with osteosarcoma's known dependency on neo-vascularization for tumor growth and metastatic colonization.
The hub genes identified through MCODE analysis point to several intertwined biological programs that collectively define the molecular character of aggressive osteosarcoma. The immune-related cluster centered on HLA-DRA and HLA-DQA1 reflects the immunological identity of the bone microenvironment. Osteoclasts, the bone-resorbing cells, are derived from the monocyte-macrophage lineage and maintain an ongoing dialogue with immune cytokines. Osteosarcoma cells appear to exploit this immune infrastructure, potentially using MHC class II upregulation to evade cytotoxic killing or to co-opt macrophage activity in ways that promote tumor survival. Dysregulation of antigen processing pathways in the NPDEGs suggests this immune disruption is an early event in osteosarcoma carcinogenesis, not a late adaptation.
Cell cycle deregulation at metastasis: The appearance of CDK1, CDK20, and CCNB1 as the highest-scoring MCODE module hub genes in the PMDEG network reinforces a well-established concept: metastatic osteosarcoma cells proliferate more aggressively than their primary tumor counterparts. The G2/M checkpoint regulated by CCNB1 and CDK1 is frequently disrupted in cancers with high metastatic potential. CDK inhibitors have entered clinical testing in multiple cancer types, and the strong network centrality of CDK1 in this osteosarcoma metastasis dataset provides disease-specific support for that therapeutic approach.
Mitochondrial translation and metabolic reprogramming: Module 2 of the PMDEG network, centered on MTIF2 (Mitochondrial Translational Initiation Factor 2) and MRPS7 (Mitochondrial Ribosomal Protein S7), implicates mitochondrial function and energy metabolism as altered in metastatic osteosarcoma. Metastatic tumor cells face significant metabolic stress during transit through the bloodstream and at colonization sites; upregulation of mitochondrial biosynthesis machinery may represent an adaptation to meet the elevated energy demands of metastatic growth.
VEGFA and angiogenesis: VEGFA (Vascular Endothelial Growth Factor A) anchors the third PMDEG module and is one of the most potent drivers of tumor angiogenesis. Its hub status in the metastatic DEG network aligns with the observation that metastatic osteosarcoma lesions require robust neo-vascularization for sustained growth at secondary sites. VEGF-targeted agents (bevacizumab and others) have been tested in osteosarcoma clinical trials, and this bioinformatics finding provides molecular rationale for continued investigation.
The authors explicitly enumerate several important limitations. First, the analysis could not account for important confounding variables present within the GEO datasets, including patient age, race, geographic region, cell lineage of origin, and tumor stage or histological subtype. These factors could introduce systematic bias into the DEG lists. Osteosarcoma encompasses multiple histological subtypes (osteoblastic, chondroblastic, fibroblastic), and pooling them without subtype stratification may obscure subtype-specific expression signatures or dilute the significance of subtype-restricted genes.
Directionality limitation: All seven genes identified as continuously deregulated across oncogenesis and metastasis were up-regulated. The authors note that continuously down-regulated genes may also exist but were not detected with sufficient consistency across all three datasets. This could reflect technical limitations of the platforms used or genuine biology, but the absence of down-regulated candidates leaves the signature incomplete. Future analyses incorporating RNA-seq datasets with greater sensitivity for low-expression transcripts may recover additional candidates.
Lack of experimental validation: The study is entirely computational. None of the identified hub genes or enriched pathways have been functionally validated in this paper using cell lines, animal models, or patient-derived samples. The SPARCL1 discrepancy with published experimental data is a concrete illustration of why computational findings require wet-lab follow-up. Validation studies using shRNA knockdown or CRISPR editing of CDK1, HLA-DRA, or VEGFA in osteosarcoma cell lines, followed by testing in xenograft mouse models, would be the natural next step.
Future directions: The authors call for prospective studies that link the expression levels of hub genes like HLA-DRA, CDK1, and VEGFA to clinical outcomes (survival, time to metastasis, chemotherapy response) in well-annotated patient cohorts. Integration with single-cell RNA-sequencing data would allow cell-type-specific expression of these genes to be resolved within the tumor microenvironment, distinguishing contributions from osteosarcoma cells, stromal cells, immune infiltrates, and endothelial cells. Drug repurposing screens targeting CDK1, VEGFA, or the TNF pathway in osteosarcoma cell lines are also suggested as immediate translational priorities.