Osteosarcoma (OS) is the most common primary malignant bone tumor, disproportionately affecting teenagers and young adults. It most frequently arises at the metaphyses of long bones, particularly around the knee, where osteoblasts are actively proliferating during skeletal growth. Despite decades of clinical research, the molecular mechanisms driving osteosarcomagenesis remain incompletely understood, and the five-year survival rate for patients with metastatic disease remains below 30%. The standard treatment, consisting of neoadjuvant chemotherapy followed by surgical resection, has changed little since the 1980s, and treatment resistance remains the dominant clinical challenge.
The role of microRNAs: MicroRNAs (miRNAs) are short, endogenous non-coding RNA molecules, typically 19-24 nucleotides in length, that regulate gene expression post-transcriptionally by binding to the 3' untranslated region (3'-UTR) of target messenger RNAs (mRNAs). This binding generally results in mRNA degradation or translational repression. A single miRNA can target dozens to hundreds of genes simultaneously, making miRNA dysregulation a high-leverage event in cancer biology. Prior studies had identified individual miRNAs with tumor-suppressive or oncogenic functions in OS, including miR-1, miR-409-3p, miR-379, miR-665, and miR-489-3p, but a systematic, network-level view of miRNA-mRNA regulatory architecture in osteosarcoma was lacking.
Study rationale: This 2017 study, published in Scientific Reports by Wang, Tang, and colleagues, leveraged publicly available gene expression data and a suite of bioinformatics tools to construct comprehensive miRNA-mRNA functional networks in osteosarcoma cell lines. Rather than focusing on individual molecules in isolation, the authors aimed to map the collective regulatory landscape, identify hub nodes in the network, and pinpoint candidate biomarkers that could be relevant for early diagnosis or targeted therapy development.
The key analytical innovation was applying an integrated pipeline using the GCBI web platform (which incorporates TargetScan and miRanda prediction databases), Cytoscape network visualization, STRING protein-protein interaction analysis, and Gene Ontology/KEGG pathway enrichment to translate raw microarray data into biologically interpretable network models.
The study drew on a single publicly deposited dataset from the Gene Expression Omnibus (GEO): accession GSE70415. This composite dataset comprises matched mRNA expression profiles (GSE70414, measured on Affymetrix Human Genome U133 Plus 2.0 Array, GPL570 platform, covering 54,675 probes) and miRNA expression profiles (GSE70367, measured on Affymetrix Multispecies miRNA-3 Array, GPL16384 platform, covering 25,533 miRNA probes). The experimental samples consisted of five established human osteosarcoma cell lines: MG63, Saos-2, HOS, NY, and Hu09. Human mesenchymal stem cells (hMSC) served as the normal comparator, representing the presumed cell of origin for osteosarcoma.
Differential expression identification: Raw expression values were processed through the GCBI platform using Median Polish normalization. Differential expression was defined by three simultaneous criteria: p-value less than 0.01 by standard statistical testing, false discovery rate (FDR) less than 0.01 using the Benjamini-Hochberg correction for multiple comparisons, and absolute fold change greater than 2 between OS cell lines and hMSC controls. These stringent thresholds were applied to both the mRNA and miRNA datasets to minimize false positives in the downstream network analysis.
Target prediction and network construction: For differentially expressed miRNAs (DEmiRNAs), predicted mRNA targets were queried from two established databases integrated within GCBI: TargetScan (seed-sequence complementarity based) and miRanda (free energy thermodynamics based). Intersection of DEmiRNA targets with the set of differentially expressed mRNAs (DEGs) defined functionally relevant miRNA-mRNA interaction pairs. Network visualization and cluster analysis were performed using Cytoscape 3.4.0, with the Molecular Complex Detection (MCODE) algorithm used to identify densely connected sub-clusters within the broader interaction network. Protein-protein interaction (PPI) networks for downstream target gene products were built using STRING 10.0.
Pathway enrichment: Gene Ontology (GO) and KEGG pathway enrichment analyses were performed on the DEG set using Fisher's exact test with FDR correction. Enrichment outputs were ranked by significance, and the top 20 GO terms and top 20 KEGG pathways were reported. Network centrality was assessed by contribution degree, a metric reflecting how many other pathways a given pathway interacts with. Only pathways with contribution degree of 10 or higher were considered as primary regulatory hubs.
Applying the differential expression criteria to the complete expression profiles, the authors identified 3,856 significantly differentially expressed genes (DEGs) across the five OS cell lines compared to hMSC controls. Of these, 1,705 were upregulated and 2,151 were downregulated in osteosarcoma, indicating that suppression of normal mesenchymal gene programs is a dominant feature of OS cell line transcriptomics. Separately, 250 significantly differentially expressed miRNAs (DEmiRNAs) were identified, with 161 showing upregulation and 89 showing downregulation in OS relative to the hMSC control.
Notable individual genes: Among the most significantly downregulated genes was Periostin (POSTN), a canonical osteoblast marker involved in cell adhesion and differentiation. The downregulation of POSTN in OS cell lines is consistent with the loss of osteoblastic differentiation capacity characteristic of osteosarcoma. However, the authors note an important caveat: prior studies using tumor biopsy specimens (rather than cell lines) have found POSTN expression to be elevated in OS tissue compared to osteochondroma, and high POSTN content in vivo correlates with tumor angiogenesis and poor prognosis. This discrepancy likely reflects differences between cell line models (which lose many in vivo signaling cues) and primary tumor specimens, as well as differences in detection platform (RNA microarray versus immunohistochemistry).
Notable miRNAs: Among the 250 DEmiRNAs, the two showing the highest absolute fold change were miR-182-5p and miR-708-5p (absolute fold change greater than 100 in both cases). Despite their extreme expression differences, neither could be found in the curated Osteosarcoma Database (which contains only 81 validated small RNA entries), highlighting the incomplete state of functional annotation for most OS-associated miRNAs. The authors note that the connection between both miRNAs and vorinostat, an FDA-approved histone deacetylase (HDAC) inhibitor, was subsequently explored in 143B and MG63 cell lines, though those results were not included in this publication.
Cross-referencing the 3,856 DEGs against the 911 curated entries in the Osteosarcoma Database revealed only approximately 7% overlap, underscoring that the vast majority of transcriptomic dysregulation in OS cell lines involves genes not yet validated by the field as OS-relevant, and that existing databases substantially underrepresent the scope of OS molecular alterations.
Functional enrichment analysis of the 3,856 DEGs identified 395 significantly enriched Gene Ontology (GO) terms and 142 significantly enriched KEGG pathways (both filtered at p < 0.01 with FDR correction). The top three enriched biological processes by GO analysis were extracellular matrix organization (GO:0030198), small molecule metabolic process (GO:0044281), and cell adhesion (GO:0007155). These three processes are tightly linked to the hallmarks of osteosarcoma biology: disrupted bone matrix remodeling, the Warburg-like metabolic reprogramming that supports rapid proliferation, and the enhanced cell-matrix interactions that facilitate local invasion and distant metastasis.
Top KEGG pathways: At the KEGG level, the three most significantly enriched pathways were: Pathways in cancer (KEGG ID: 5200, degree 27), PI3K-Akt signaling pathway (KEGG ID: 4151), and Metabolic pathways (KEGG ID: 1100). The predominance of metabolic pathway enrichment supports the emerging understanding that chemoresistance in OS is partly mediated by metabolic abnormalities, consistent with studies showing that miR-221, miR-101, miR-22, and miR-155 participate in cisplatin and doxorubicin resistance via metabolic and autophagy mechanisms.
Pathway network centrality: Among all enriched KEGG pathways, the MAPK signaling pathway (degree 34), Pathways in cancer (degree 27), and Cell cycle (degree 24) showed the highest contribution degrees in the co-enrichment network, meaning they interact with the greatest number of other significantly enriched pathways. This positions MAPK, general cancer pathway cascades, and cell cycle regulation as the most central regulatory hubs mediating downstream pathway perturbation in OS. The p53 signaling pathway (degree 17), Wnt signaling (degree 16), TGF-beta signaling (degree 13), and Focal adhesion (degree 12) also exhibited high centrality.
The PI3K-Akt pathway, second in KEGG enrichment significance, is particularly relevant clinically because its activation suppresses FOXO transcription factors via phosphorylation, promoting cell survival. Crosstalk between PI3K-Akt and MAPK pathways, as well as NF-kB signaling, creates a reinforcing network of pro-survival signaling. Activation of MAPK signaling through elevated EGFR phosphorylation and MMP-9 levels, mediated in part by loss of the tumor-suppressive miR-143, represents one specific mechanistic axis captured in this enrichment analysis.
To construct the miRNA-mRNA regulatory network, the authors queried predicted targets for the 250 DEmiRNAs from the integrated TargetScan/miRanda databases within GCBI, identifying 29,227 genes deposited across the two target pools. Of these, 388 were substantially involved in GO-enriched functional categories, and 608 overlapped with the set of differentially expressed mRNAs regardless of binding pair specificity. To focus on functionally meaningful interactions, only miRNAs with an interaction degree of 10 or higher (meaning they regulate at least 10 mRNA targets) were considered significant for network visualization. Forty DEmiRNAs meeting this threshold were selected for detailed network mapping.
Network results: The resulting miRNA-mRNA interaction networks comprised 1,181 regulatory linkages in total. Within the network of upregulated miRNAs, 238 downstream target genes were found to be in a repressed state (consistent with miRNA-mediated suppression), while in the network of downregulated miRNAs, 181 targets were found to be in an activated state (consistent with derepression when their miRNA regulators are lost). This bidirectional regulatory architecture reflects the simultaneous gain of oncomiR activity and loss of tumor-suppressive miRNA function in osteosarcoma.
Hub miRNAs: Two miRNAs emerged as the dominant regulatory hubs based on interaction degree, meaning they showed the most significant collective impact on gene transcription across the network. miR-93-5p, which is overexpressed in OS, was identified as the primary oncomiR hub, with the highest number of downstream suppressive interactions. miR-29b-3p, which is underexpressed in OS, emerged as the primary tumor-suppressive hub. Previous work had confirmed that miR-29b-3p exerts tumor-suppressive functions in OS with tumor-specific subcellular localization patterns. Notably, the two miRNAs with the largest absolute fold changes, miR-182-5p and miR-708-5p, showed only moderate contribution degrees in the interaction network, reinforcing the point that the degree of expression change does not necessarily correlate with the breadth of regulatory impact.
The authors used Cytoscape 3.4.0 to visualize both the upregulated miRNA regulatory network and the downregulated miRNA regulatory network separately, representing miRNAs as triangles or circles (colored red) and their target mRNAs as blue squares. This visual separation of oncomiR and tumor-suppressive miRNA networks provides a clearer mechanistic framework than a single aggregate network.
Beyond the miRNA-mRNA regulatory networks, the authors performed two complementary analyses to further characterize functional organization within the OS transcriptome. First, they constructed a gene co-expression network using 698 overlapping genes derived from the intersection of GO-enriched and KEGG-enriched DEGs. In a co-expression network, edges connect genes whose expression levels change together across samples, suggesting either direct regulatory relationships or membership in the same functional module, even without a defined upstream regulatory mechanism.
CKMT2 as the top co-expression hub: Application of the MCODE (Molecular Complex Detection) algorithm to the co-expression network identified the most densely connected sub-cluster, comprising 19 closely associated genes, all of which were intimately connected to creatine kinase, mitochondrial 2 (CKMT2), also known as SMTCK. CKMT2 is an enzyme that plays an indispensable role in maintaining rational energy metabolism through the creatine phosphate shuttle, transferring high-energy phosphate from mitochondria to cytosolic sites of energy demand. Its emergence as the central hub of the most tightly co-expressed gene module in OS suggests a critical role for mitochondrial energy metabolism dysregulation in OS biology, a finding the authors indicate was subsequently validated by a collaborating group (though unpublished at the time of this paper).
Protein-protein interaction network: For the protein-protein interaction analysis, the authors focused on 35 DEmiRNAs with fold change of 10 or higher (a stricter subset than the 40 used for miRNA-mRNA network mapping). The mRNA targets of these high-magnitude DEmiRNAs were submitted to STRING 10.0 for PPI network construction, yielding 43 genes for Cytoscape-based visualization. Within this PPI network, several classical tumor suppressor and oncogenic proteins emerged as interaction nodes, including TP53, EGFR, MMP2, FOXO1, BMP family members, and members of the collagen (COL) and integrin (ITG) families. The convergence of miRNA regulatory influence on TP53, EGFR, and MMP family members provides a mechanistic link between the non-coding RNA dysregulation and the canonical oncogenic pathways that drive OS invasion, angiogenesis, and metastasis.
The two hub miRNAs identified in this network analysis, miR-29b-3p and miR-93-5p, have individually been the subject of prior experimental studies that help contextualize the network-level findings. miR-29b-3p is consistently found downregulated across multiple OS studies and exerts tumor-suppressive effects through a range of downstream targets. Its loss derepresses extracellular matrix components including collagen family members, contributing to the matrix remodeling that facilitates local invasion. Prior work had validated miR-29b-3p tumor suppression in OS with tumor-specific subcellular localization, providing a mechanistic basis for why its loss is selectively advantageous during OS development.
miR-93-5p as an oncomiR: miR-93-5p, a member of the miR-17-92 cluster family, promotes proliferation and survival in multiple cancer types. In OS, its overexpression contributes to the activation of PI3K-Akt signaling through suppression of negative regulators within that pathway. The PI3K-Akt axis in turn phosphorylates and inactivates FOXO transcription factors, removing a key brake on cell proliferation and survival. miR-93-5p also participates in regulating EGFR pathway activity, linking it to the MAPK signaling enrichment identified in the pathway analysis. The convergent impact of miR-93-5p on PI3K-Akt-MAPK signaling makes it a plausible therapeutic target, though no inhibitors of this specific miRNA have reached clinical evaluation in OS.
Chemoresistance implications: Several of the miRNAs identified in this dataset, including miR-221, miR-101, miR-22, and miR-155, have been independently validated as mediators of chemoresistance in OS. miR-221 activates PI3K-Akt to promote cisplatin resistance. miR-101 promotes chemosensitivity by blocking autophagy. miR-22 inhibits autophagy through suppression of HMGB1. miR-155 promotes drug resistance through autophagy induction. The simultaneous dysregulation of multiple resistance-associated miRNAs in OS cell lines supports the hypothesis that miRNA network-level changes, rather than individual miRNA alterations, drive the multidrug resistance phenotype characteristic of relapsed OS.
EGFR and MAPK connectivity: Within the PPI network, miR-143 loss was specifically linked to elevated EGFR phosphorylation and MMP-9 upregulation, mechanistically connecting a single miRNA depletion event to two separate established drivers of OS invasion: receptor tyrosine kinase hyperactivation and matrix metalloproteinase-mediated extracellular matrix degradation. These two processes are foundational to the high metastatic potential that defines clinical OS behavior.
Cell line versus primary tumor limitations: The most consequential limitation of this study is that all expression data derive from cell lines rather than primary patient tumor specimens. As illustrated by the POSTN discrepancy, cell lines can show expression patterns that are directionally opposite to those observed in actual tumor tissue. Cell lines lack the complex microenvironmental signals, stromal components, immune infiltration, and extracellular matrix architecture of in vivo tumors. Findings validated in cell lines require independent confirmation in patient-derived xenografts and primary tissue cohorts before they can be considered biologically relevant to human OS disease.
Single dataset dependency: The entire analysis rests on a single GEO dataset (GSE70415) comprising only five OS cell lines and one control cell type (hMSC). With such a small sample size (n=5 per group), the statistical power to detect subtle differential expression is limited, and the results may be heavily influenced by the specific biology of these particular cell lines rather than representing general OS transcriptomics. The Osteosarcoma Database overlap rate of only 7% further highlights how much of the identified transcriptomic signature remains unvalidated.
Prediction-based miRNA targets: The miRNA-mRNA interactions in this network are almost entirely based on computational target prediction (TargetScan, miRanda) rather than experimental validation (e.g., luciferase reporter assays, Argonaute immunoprecipitation, or CLASH sequencing). Target prediction algorithms generate substantial false-positive rates, meaning many of the 1,181 linkages in the network may not represent genuine regulatory relationships. Experimental validation of even a subset of the identified hub miRNA-mRNA pairs would substantially strengthen the biological conclusions.
Future directions: The authors identify two priority extensions: (1) incorporating additional non-coding RNA classes, particularly long non-coding RNAs (lncRNAs) and circular RNAs (circRNAs), which participate in competitive endogenous RNA (ceRNA) mechanisms and add another regulatory layer above the miRNA-mRNA axis studied here, and (2) integrating small molecule inhibitor and drug data to connect the network to actionable therapeutic targets. The authors also note that the CKMT2 hub finding and the connection between miR-182-5p/miR-708-5p and the HDAC inhibitor vorinostat were being actively pursued as follow-up projects at the time of publication. More broadly, the bioinformatics pipeline demonstrated here, combining GCBI, Cytoscape, and STRING, represents a replicable approach that can be applied to larger patient cohort datasets as they become publicly available.