Osteosarcoma (OS) is the most common non-haematological primary malignant bone tumor, occurring predominantly in the metaphyseal regions of long bones in adolescents and young adults, as well as in patients over 40. Despite significant improvement in survival following the introduction of neoadjuvant chemotherapy, the overall prognosis remains poor, particularly for patients who develop metastatic or chemotherapy-resistant disease.
OS is characterized by complex and highly variable karyotypes with numerous genomic aberrations. This genomic complexity means that individual gene expression studies often produce inconsistent results due to small sample sizes, different platforms, and heterogeneous patient populations. Single-study findings are thus difficult to replicate and may reflect noise rather than true biological signals.
Meta-analysis, the systematic integration of gene expression data from multiple independent studies, offers a solution by increasing statistical power, allowing heterogeneity assessment, and producing more robust and reproducible findings. This study performed the first meta-analysis of OS microarray gene expression data to identify consistently differentially expressed genes (DEGs) and the biological processes they represent.
The authors searched PubMed and the Gene Expression Omnibus (GEO) database using the terms osteosarcoma, gene expression, microarray, and genetics. After applying inclusion criteria (only original experiments comparing OS to normal control tissues, human studies only), 8 datasets were retained. In total, the meta-analysis included 240 OS samples and 35 normal control samples.
The datasets were profiled on multiple platforms: two used Affymetrix HG-U133A arrays (GPL96), two used Affymetrix HuGene-1_0-st arrays (GPL6244), two used Illumina human-6 v2.0 beadchips, one used Illumina HumanHT-12 V3.0, and one used Affymetrix HG-U133_Plus_2. Sample sources included 4 studies of in vivo bone tissues, 3 studies of OS cell lines in vitro, and 1 combined study. GEO accession numbers ranged from GSE11414 to GSE42352.
Data from these diverse sources were harmonized using a Z-score normalization approach, transforming raw probe intensities to z-scores using the formula Z = (xi - mean) / SD within each experiment. This global normalization minimized cross-platform technical variability while preserving within-dataset biological signal.
After Z-score transformation, all datasets were merged and the Significance Analysis of Microarrays (SAM) method was applied to identify differentially expressed genes. SAM uses gene-specific t-statistics with a relative difference score (D value) defined as the average expression change divided by the standard deviation of measurements for that gene. Genes were selected as significantly differentially expressed if they showed at least a 2-fold change with a false discovery rate (FDR) below 0.05.
Functional annotation was performed using GENECODIS for Gene Ontology (GO) enrichment analysis across three categories: biological process, molecular function, and cellular component. Genes below a nominal significance threshold of p less than 0.01 were tested against the background set of all annotated genes. KEGG pathway enrichment was also performed using a hypergeometric test with p less than 0.05 as the selection criterion.
Protein-Protein Interaction (PPI) network construction used the BioGRID database, with visualization in Cytoscape. The PPI analysis focused specifically on the top 10 most significantly up-regulated and down-regulated DEGs, producing a network with 129 nodes and 182 edges to identify hub proteins with the highest degree of connectivity.
The meta-analysis identified 979 DEGs across the 8 studies: 472 up-regulated and 507 down-regulated in OS compared to normal controls. The most significantly up-regulated gene was CPE (carboxypeptidase E), with a p-value of 5.08 x 10^-15 and a fold change of 3.856. CPE is involved in the biosynthesis of peptide hormones and neurotransmitters including insulin, and has previously been correlated with tumor growth and metastasis in pheochromocytomas and other cancers, though its role specifically in OS had not been previously reported.
Other notable up-regulated genes in the top 10 included HEY1 (p = 6.72 x 10^-15, fold change 2.71), a NOTCH1 target gene previously found elevated in murine and canine OS; FXYD6 (p = 2.40 x 10^-11, fold change 2.39), previously enriched in OS cell lines; and EFNA1 (ephrin-A1, p = 1.23 x 10^-10, fold change 1.95), implicated in MAPK activation via the EphA2 receptor in OS. The most significantly down-regulated gene was NPR3 (natriuretic peptide receptor 3, p = 1.86 x 10^-48, fold change -2.35), followed by GAS6 (p = 2.5 x 10^-31, fold change -2.69) and RGS4 (p = 8.3 x 10^-27, fold change -3.26).
GAS6 (growth arrest-specific 6) is particularly biologically relevant: in OS cell lines, recombinant human GAS6 activates the Axl receptor tyrosine kinase to protect tumor cells from apoptosis induced by serum starvation, and promotes tumor cell migration and invasion in vitro, suggesting that its down-regulation may reflect a compensatory suppressive mechanism or represent a loss of function with complex consequences.
GO enrichment analysis identified significant enrichment in protein binding (GO: 0005515, p = 3.83 x 10^-60) and calcium ion binding (p = 3.79 x 10^-13) for molecular functions; cell adhesion (p = 2.26 x 10^-19) and negative regulation of apoptotic process (p = 3.24 x 10^-15) for biological processes; and cytoplasm (p = 9.18 x 10^-63) and extracellular region (p = 2.28 x 10^-47) for cellular components.
KEGG pathway analysis revealed that the most significantly enriched pathway was Focal adhesion (hsa04510, p = 5.70 x 10^-15, 34 genes). This is followed by ECM-receptor interaction (p = 1.27 x 10^-13, 22 genes) and Cell cycle (p = 4.53 x 10^-11, 23 genes). Other enriched pathways included cytokine-cytokine receptor interaction, complement and coagulation cascades, regulation of actin cytoskeleton, and osteoclast differentiation.
The focal adhesion pathway finding is particularly notable because focal adhesions are associated with cell migration dynamics, which appears paradoxical given the migratory phenotype expected in metastatic OS. However, prior work demonstrated that knockdown of paxillin in highly metastatic OS sub-lines M112 and 132 inhibits cell migration, resolving this apparent contradiction and confirming that focal adhesion signaling genuinely contributes to OS invasiveness.
The PPI network constructed from the top 10 up-regulated and down-regulated DEGs in BioGRID comprised 129 nodes and 182 edges. The three hub proteins with the highest connectivity (degree) were PTBP2 (polypyrimidine tract binding protein 2, degree = 33), RGS4 (regulator of G-protein signaling 4, degree = 15), and FXYD6 (FXYD domain containing ion transport regulator 6, degree = 13).
PTBP2 is a member of the polypyrimidine tract binding protein family, which regulates post-transcriptional events including alternative splicing. It is expressed in the nervous system, neural retina, spinal cord, and intermediate mesoderm, and represses adult-specific splicing to regulate neuronal precursor generation in the embryonic brain. Its function in OS had not been described, and its high connectivity in the PPI network suggests it may be a previously unrecognized OS driver gene.
RGS4 (down-regulated, fold change -3.26) is a regulator of G-protein signaling and ranked as the most significantly down-regulated gene in the PPI hub set. FXYD6 (up-regulated, fold change 2.39), an ion transport regulator, was previously identified as enriched in OS cell lines by directional tag PCR subtractive hybridization. These hub proteins represent high-priority candidates for mechanistic follow-up studies.
The fundamental challenge of this meta-analysis is the heterogeneity of the included datasets. Clinical samples may differ in disease activity, stage, gender distribution, and treatment history. Different microarray platforms use different probe designs and normalization algorithms, and while Z-score normalization reduces this variability, it cannot fully eliminate platform-specific biases.
A significant confounding factor is that data from 4 studies were derived from in vivo bone tissue, 3 from OS cell lines, and 1 from both sources. Gene expression in cell lines may diverge substantially from primary tumor biology due to culture adaptation, loss of stromal context, and passage-related changes. These tissue-type differences were not formally modeled as covariates in the analysis.
The meta-analysis detected consistent signals across studies despite this heterogeneity, which argues for robustness of the top DEGs. However, the authors acknowledge that further experimental research is needed to confirm findings. The small number of control samples (35 total) relative to OS samples (240) may also limit the precision of fold change estimates for less robustly regulated genes.
The identification of CPE, HEY1, FXYD6, and EFNA1 as consistently up-regulated genes in OS provides a focused list of candidates for functional validation. Notch signaling via HEY1 is a particularly actionable target, as gamma-secretase inhibitors and anti-Notch biologics are already in clinical development for various cancers and could be evaluated in OS models.
PTBP2's high connectivity in the PPI hub network and its known role in splicing regulation suggest that OS may exploit alternative splicing programs, analogous to mechanisms described in brain cancer. Systematic splicing analysis of OS RNA-seq data using PTBP2 as an anchor could reveal novel OS-specific splice isoforms relevant to diagnosis or therapy.
Future meta-analyses should incorporate RNA-sequencing datasets as they become available, which would provide allelic resolution and reveal non-coding RNA dysregulation beyond what microarray platforms can detect. Integration of DNA methylation, copy number variation, and proteomic data from the same cohorts would enable a truly multi-omics view of the OS molecular landscape.