Osteosarcoma is the most common primary malignant bone tumor in children and adolescents, with a worldwide average incidence of 3.1 per million population overall and 4.4 per million in those under 25 years of age. The tumor shows a bimodal age distribution, peaking in adolescence and again in older adults. The standard treatment protocol combines neoadjuvant chemotherapy (NACT) before surgery with adjuvant chemotherapy afterward. Prior to the introduction of multiagent chemotherapy in the 1970s, fewer than 20% of patients survived five years; that figure climbed to 60-70% for localized disease once systemic treatment was adopted.
The chemoresistance crisis: Despite this progress, five-year survival for patients who respond poorly to initial chemotherapy remains below 25%. The 40% of patients who are primary non-responders, or who relapse after achieving an initial response, face a dismal prognosis. Current clinical tools cannot reliably identify in advance which patients will be refractory to doxorubicin- and cisplatin-based regimens, and the molecular mechanisms driving chemoresistance in osteosarcoma are only partially understood.
The ABCB1 precedent: The best-characterized chemoresistance gene in osteosarcoma encodes multidrug resistance protein 1 (MDR1), the product of the ABCB1 gene. MDR1 functions as an ATP-dependent drug efflux pump and has been validated across multiple osteosarcoma studies as a correlate of doxorubicin resistance and adverse outcome. The existence of ABCB1 as a validated target motivates the search for additional, potentially more druggable molecular drivers of chemoresistance that have been missed by earlier candidate-gene approaches.
This 2018 study published in BioMed Research International applies a systematic bioinformatics pipeline to two publicly available Gene Expression Omnibus (GEO) datasets to identify differentially expressed genes (DEGs) and microRNAs (DEMs) between osteosarcoma patients with poor versus good chemotherapy responses. The goal is to nominate novel hub genes and miRNA targets that could serve as therapeutic targets or prognostic biomarkers.
The study draws on two datasets retrieved from the National Center for Biotechnology Information Gene Expression Omnibus (NCBI GEO) database. GSE87437 is a gene expression microarray dataset based on the GPL570 platform (Affymetrix HG-U133 Plus 2.0 Array), containing samples from 10 poor-response and 11 good-response osteosarcoma patients. GSE30934 is a miRNA expression array based on the GPL10312 platform (3D-Gene Human miRNA Oligo chip v12-1.00), containing 8 poor-response and 16 good-response samples. These datasets were submitted by independent research groups, providing some degree of methodological independence.
Differential expression analysis: GEO2R, an R-based web application provided by NCBI, was used to identify statistically significant expression differences between the two groups in each dataset. The cutoff criteria applied were P less than 0.05 and an absolute log2 fold change of 1.0 or greater. This two-threshold approach requires both statistical significance and biological magnitude, reducing false positives from the large number of probes tested simultaneously. A hierarchical clustering analysis of the top 100 DEGs was visualized using the Morpheus web tool to confirm that differential expression separated good-response from poor-response samples.
Functional enrichment analysis: All 668 DEGs were uploaded to the Database for Annotation, Visualization, and Integrated Discovery (DAVID v6.8) for Gene Ontology (GO) and Kyoto Encyclopedia of Genes and Genomes (KEGG) pathway enrichment analysis. GO analysis covered three domains: biological process (BP), molecular function (MF), and cellular component (CC). KEGG analysis linked the DEG set to known signaling and metabolic pathways at a significance threshold of P less than 0.05.
Protein-protein interaction network and hub gene selection: DEGs were submitted to the STRING database (covering 9,643,763 proteins from 2,031 organisms) with a minimum interaction confidence score of 0.4. The resulting PPI network was visualized in Cytoscape and analyzed using the Molecular Complex Detection (MCODE) plugin to identify tightly connected subnetworks or "modules." MCODE parameters were: degree cutoff = 2, node density cutoff = 0.1, node score cutoff = 0.2, k-core = 2, max depth = 100. Hub genes were exported from the top modules, and their prognostic significance was assessed using GSE21257, a separate GEO dataset containing osteosarcoma patient survival data.
Applying the dual thresholds of P less than 0.05 and |log2FC| of 1.0 or greater to GSE87437 yielded a total of 668 differentially expressed genes between poor-chemotherapy-response and good-chemotherapy-response osteosarcoma samples. Of these, 422 were upregulated in poor responders, meaning their expression was higher in chemoresistant tumors, and 246 were downregulated. The asymmetry between upregulated and downregulated genes suggests that chemoresistance is associated with a net increase in transcriptional activity across a broad range of biological processes rather than a uniform suppression.
Hierarchical clustering: When the top 100 DEGs (50 most upregulated and 50 most downregulated) were used for unsupervised hierarchical clustering, the resulting heat map showed clear separation between good-response and poor-response samples. This visual confirmation that gene expression differences are large enough to cluster samples by phenotype is an important internal validation step, demonstrating that the identified DEGs carry genuine biological signal rather than statistical noise driven by the relatively small sample sizes.
The 668 DEGs represent a shortlist from a much larger set of probed transcripts on the Affymetrix HG-U133 Plus 2.0 platform, which interrogates over 54,000 probe sets. Identifying 668 DEGs from a 21-sample dataset is consistent with the expected yield for a fold-change-plus-p-value filtering approach applied to a biologically meaningful phenotypic contrast. The downstream GO, KEGG, and PPI analyses are designed to distill this list further into actionable biological insights.
Gene Ontology enrichment analysis categorized the 668 DEGs across three annotation domains. In the biological process domain, the most significantly enriched terms were positive regulation of transcription (DNA-templated), positive regulation of sequence-specific DNA binding transcription factor activity, nitric oxide-mediated signal transduction, positive regulation of transcription from RNA polymerase II promoter, and regulation of phosphatidylinositol 3-kinase (PI3K) signaling. The preponderance of transcriptional regulation terms suggests that chemoresistance in osteosarcoma is, at least in part, driven by rewiring of gene expression programs rather than by mutations in individual drug-metabolizing enzymes.
Molecular function and cellular component: In the molecular function domain, enriched terms included estrogen response element binding, Rac guanyl-nucleotide exchange factor activity, calcium ion binding, zinc ion binding, and phosphatidylinositol-4,5-bisphosphate 3-kinase (PIP2 kinase) activity. Cellular component enrichment highlighted proteinaceous extracellular matrix, cell surface, P granule, integral component of plasma membrane, and endocytic vesicle membrane. The combined picture implicates altered extracellular matrix composition and membrane-associated signaling as components of the chemoresistant phenotype.
KEGG pathway analysis: KEGG enrichment identified five principal pathways: tryptophan metabolism, oxytocin signaling pathway, glyoxylate and dicarboxylate metabolism, cyclic AMP (cAMP) signaling pathway, and dopaminergic synapse. The tryptophan metabolism pathway is of particular interest because tryptophan catabolism via the kynurenine pathway is a well-characterized immune-evasion mechanism in tumors, and its upregulation in chemoresistant osteosarcoma samples could link treatment resistance to simultaneous immune escape. The cAMP signaling pathway has established roles in modulating drug sensitivity across cancer types through downstream regulation of pro-survival kinases and transcription factors.
Module-specific KEGG analysis of the three PPI subnetworks revealed additional pathway associations: ribosome biogenesis in eukaryotes (module 1), calcium signaling pathway and arachidonic acid metabolism (module 2), and proteoglycans in cancer and linoleic acid metabolism (module 3). These module-level pathways provide more targeted mechanistic hypotheses for each functional cluster of interacting proteins.
STRING-based PPI network construction using the 668 DEGs produced a network of 432 nodes connected by 428 edges. The moderate edge density relative to node count indicates a network where most proteins interact with only a small number of partners, typical of biological interaction networks that follow a power-law degree distribution. After importing this network into Cytoscape and applying MCODE module detection, three significant subnetworks were identified, each representing a functionally cohesive cluster of proteins with elevated internal connectivity.
The nine hub genes: From these modules, nine hub genes were extracted: ZNRD1 (zinc ribbon domain containing 1), MYH7B (myosin heavy chain 7B), GPR68 (G protein-coupled receptor 68), CAT (catalase), FUT3 (fucosyltransferase 3, Lewis blood group), IMPG2 (interphotoreceptor matrix proteoglycan 2), GPR180 (G protein-coupled receptor 180), ANPEP (alanyl aminopeptidase, membrane), and CDK1 (cyclin dependent kinase 1). Seven of the nine, including ZNRD1, MYH7B, GPR68, FUT3, IMPG2, ANPEP, and CDK1, were upregulated in chemoresistant samples; two, CAT and GPR180, were downregulated.
Functional significance of the hub genes: CDK1 is a master regulator of cell-cycle progression whose upregulation in chemoresistant tumors is consistent with accelerated proliferation and reduced sensitivity to cytotoxic agents. ANPEP encodes a membrane aminopeptidase known to be silenced in prostate cancer, and its upregulation in chemoresistant osteosarcoma may reflect an altered protease landscape facilitating invasion and drug efflux. GPR68 is a pH-sensing G protein-coupled receptor linked to tumor microenvironment acidification, a property associated with multidrug resistance because acidic conditions reduce intracellular drug accumulation. CAT encodes catalase, the primary enzymatic scavenger of hydrogen peroxide; its downregulation in poor responders is consistent with prior work showing that reduced catalase activity promotes reactive oxygen species accumulation and, paradoxically, activates pro-survival signaling cascades that blunt drug-induced apoptosis.
To evaluate whether the nine hub genes carry independent prognostic information, the authors performed Kaplan-Meier survival analysis using GSE21257, a separate GEO dataset that includes osteosarcoma patient survival data. Note that survival data for GPR180 and CDK1 were not available within GSE21257, so these two hub genes could not be assessed for prognostic impact using this dataset. For the remaining seven genes with available data, the most striking finding involved FUT3.
The FUT3 paradox: FUT3 (fucosyltransferase 3, Lewis blood group) encodes an enzyme that catalyzes the addition of fucose sugar residues to glycoprotein substrates, producing sialyl Lewis X and related carbohydrate antigens on cell surface proteins. FUT3 expression was significantly upregulated in chemoresistant osteosarcoma samples in the GSE87437 analysis, classifying it as a driver or marker of poor treatment response. Yet when survival analysis was performed in GSE21257, osteosarcoma patients with high FUT3 mRNA expression showed better overall survival compared to those with low expression.
Interpreting the paradox: This seemingly contradictory finding may reflect the complexity of sialyl Lewis X biology in cancer. On one hand, fucosylation promotes cancer stem cell invasiveness, metastasis, and immune evasion. On the other hand, sialyl Lewis X antigens serve as ligands for natural killer (NK) cell lectin-like receptors, and multiple studies have shown that NK cells exhibit enhanced cytotoxicity against tumor cells expressing high sialyl Lewis X levels. If high FUT3 expression simultaneously marks chemoresistant cells while also rendering them more visible to NK-mediated immune surveillance, the net survival effect could appear paradoxically favorable. This finding highlights the need for careful, context-specific interpretation of bioinformatics-derived biomarkers before clinical translation.
Parallel analysis of the miRNA expression array GSE30934 using the same GEO2R pipeline and cutoff criteria (P less than 0.05, |log2FC| greater than or equal to 1.0) identified five differentially expressed microRNAs (DEMs) between poor-response and good-response osteosarcoma patients. MicroRNAs are short non-coding RNAs of approximately 18-25 nucleotides that suppress gene expression post-transcriptionally by binding to complementary sequences in the 3' untranslated regions (3' UTR) of target mRNAs, either blocking translation or triggering mRNA degradation.
The five DEMs: hsa-miR-543 (log2FC -3.43, P = 0.00192), hsa-miR-409-5p (log2FC -2.71, P = 0.00729), hsa-miR-154 (log2FC -2.62, P = 0.03838), hsa-miR-518f (log2FC +1.45, P = 0.02332), and ebv-miR-BART1-3p (log2FC -1.36, P = 0.04733). Four of the five miRNAs were downregulated in poor responders, with hsa-miR-543 showing the largest magnitude of change (log2FC -3.43). hsa-miR-518f was the only upregulated miRNA, suggesting it may be playing a pro-resistance role. The presence of ebv-miR-BART1-3p, a miRNA encoded by Epstein-Barr virus, raises the question of whether viral co-infection contributes to the chemoresistant phenotype in some osteosarcoma patients.
miRNA-DEG network construction: To link DEMs with DEGs, the authors used miRWalk 1.0, an integrated platform querying ten miRNA target prediction databases simultaneously (DIANAmT, miRanda, miRDB, miRWalk, RNAhybrid, PICTAR4, PICTAR5, PITA, RNA22, and Targetscan). The strength of predicted miRNA-target interactions was color-coded by the number of databases supporting each pair, with red indicating the strongest cross-database consensus. Cross-referencing the predicted targets with the nine hub genes identified two high-confidence regulatory pairs: hsa-miR-543 targeting ZNRD1, and hsa-miR-518f targeting CAT. Both pairs showed concordant expression trends, meaning the miRNA and its putative target gene moved in opposite directions between good and poor responders, consistent with genuine regulatory relationships.
Small sample sizes: The most significant limitation of this study is the small number of patient samples in each GEO dataset. GSE87437 contains only 21 patients and GSE30934 contains only 24 patients. At these sample sizes, the statistical power to reliably detect all truly differentially expressed genes is limited, and false-positive findings are likely mixed into the 668 DEGs and 5 DEMs. The reliance on two existing publicly deposited datasets also means the authors had no control over patient selection, treatment regimens, chemotherapy response criteria, or biopsy timing, all of which can introduce heterogeneity that obscures biologically meaningful signals.
Computational prediction vs. experimental evidence: The miRNA-DEG regulatory pairs (hsa-miR-543/ZNRD1 and hsa-miR-518f/CAT) were identified through computational prediction databases rather than direct experimental validation. While the concordance of expression trends is encouraging, prediction algorithms for miRNA targets have variable precision, and many computationally predicted miRNA-target interactions fail to validate in functional assays. Luciferase reporter assays, miRNA overexpression and knockdown experiments in osteosarcoma cell lines, and xenograft models are needed to confirm that these miRNAs actually regulate their predicted targets and modulate chemosensitivity.
Lack of mechanistic depth: The study does not investigate the mechanisms by which individual hub genes drive chemoresistance. For most of the nine hub genes, the available literature on their osteosarcoma-specific biology was sparse at the time of publication. For example, MYH7B was primarily characterized in cardiac biology, and IMPG2 was known as a retinal proteoglycan associated with retinitis pigmentosa, with no prior connection to chemoresistance. Whether these genes are causal drivers of resistance or correlative bystanders co-regulated with the resistance phenotype cannot be determined from expression data alone.
Future directions: The authors propose that the identified hub genes and miRNA pairs represent a prioritized set of candidates for functional experimentation. Key next steps include in vitro validation using osteosarcoma cell lines with established chemosensitive and chemoresistant phenotypes, in vivo validation in mouse xenograft models, and ultimately clinical validation in larger, prospectively collected patient cohorts. The miRNA-target pairs are particularly attractive therapeutic candidates because miRNA mimics and inhibitors are clinically tractable molecules currently in development for multiple cancer types.