IJMM International Journal of Molecular Medicine 1107-3756 1791-244X D.A. Spandidos 10.3892/ijmm.2019.4108 ijmm-43-04-1723 Articles Identification of biomarkers associated with the recurrence of osteosarcoma using ceRNA regulatory network analysis ZhangShanyong1 DingLei2 LiXin3 FanHongwu4 Departments of Spinal Surgery Rehabilitation Joint Surgery, The Second Hospital of Jilin University, Changchun, Jilin 130041 Department of Orthopedics, China Japan Union Hospital of Jilin University, Changchun, Jilin 130021, P.R. China Correspondence to: Dr Hongwu Fan, Department of Orthopedics, China Japan Union Hospital of Jilin University, 126 Xiantai Street, Changchun, Jilin 130021, P.R. China, E-mail: hongwwwwfff@outlook.com 04 2019 25 02 2019 43 4 1723 1733 27 03 2018 07 12 2018 Copyright: © Zhang et al. 2019 This is an open access article distributed under the terms of the Creative Commons Attribution-NonCommercial-NoDerivs License, which permits use and distribution in any medium, provided the original work is properly cited, the use is non-commercial and no modifications or adaptations are made.

The aim of the present study was to identify the important mRNAs, micro (mi)RNAs and long non-coding (lnc)RNAs that are associated with osteosarcoma recurrence. The GSE3905 dataset, which contains two sub-datasets (GSE39040 and GSE39055), was downloaded from the Gene Expression Omnibus (GEO). Prognosis-associated RNAs were identified by performing Cox regression univariate analysis and were subsequently used to construct a competing endogenous (ce)RNA regulatory network for Gene Set Enrichment Analysis (GSEA). Kaplan-Meier survival analysis was used to determine the associations between expression levels and survival prognosis. In addition, another independent miRNA profile, GSE79181, was downloaded from GEO for validation. Among the differentially expressed RNAs, 417 RNAs (5 lncRNAs, 19 miRNAs, and 393 mRNAs) were observed to be associated with prognosis. The GSEA for the ceRNA regulatory network revealed that ‘Mitogen-activated protein kinase (MAPK) signaling pathway’, ‘Chemokine signaling pathway’ and ‘Spliceosome’ were markedly associated with osteosarcoma. In addition, three lncRNAs [long intergenic non-protein coding RNA 28 (LINC00028), LINC00323, and small nucleolar RNA host gene 1 (SNHG1)] and two miRNAs (hsa-miR-124 and hsa-miR-7) regulating three mRNAs [Ras-related protein Rap-1b (RAP1B), activating transcription factor 2 (ATF2) and protein phosphatase Mg2+/Mn2+ dependent 1B (PPM1B)] participated in the MAPK signaling pathway. The Kaplan-Meier survival analysis also demonstrated that samples with lower expression levels of LINC00323 and SNHG1 had better prognosis, and samples with increased expression levels of LINC00028, hsa-miR-124 and hsa-miR-7 had better prognosis. Overexpression of RAP1B, ATF2 and PPM1B was positively associated with osteosarcoma recurrence. The roles of hsa-miR-124 and hsa-miR-7 in osteo-sarcoma recurrence were also validated using GSE79181. Thus, in conclusions, the three lncRNAs (LINC00028, LINC00323 and SNHG1), two miRNAs (hsa-miR-124 and hsa-miR-7) and three mRNAs (RAP1B, ATF2, and PPM1B) were associated with osteosarcoma recurrence.

 osteosarcoma recurrence competing endogenous RNA network Introduction

Osteosarcoma is the most common primary malignant bone tumor and often exhibits locally invasive growth. Osteosarcoma is commonly characterized by an appendicular primary tumor with a high rate of pulmonary metastases (1). In addition, osteosarcoma is a genetically unstable and highly malignant mesenchymal tumor of the bone characterized by structural chromosomal alterations (2,3).

Improved understanding of the molecular mechanisms underlying osteosarcoma would support the development of effective and targeted therapies. Many previous studies have reported molecular mechanisms underlying osteosarcoma that were involved in the progression and metastasis of osteosarcoma (4-8). Inhibiting Wnt/β-catenin signaling and targeting c-Met, a Wnt-regulated proto-oncogene, have shown to be useful in the treatment of osteosarcoma (4,5). The Notch pathway has a vital role in regulating tumor angiogenesis and vasculogenesis in osteosarcoma (6). Wilms' tumor gene 1 has been reported to act as a tumor promoter in osteosarcoma by regulating cell cycle and apoptotic machinery (7). Recently, COP9 signalosome subunit 3 (COPS3), another oncogene overexpressed in osteosarcoma, was revealed to be involved in osteosarcoma metastasis by interacting with Beclin1 and Raf-1 through autophagy (8). In addition, some studies have reported the roles of long non-coding RNAs (lncRNAs) in the progression, metastasis, and prognosis of osteosarcoma (9-16). A previous study also demonstrated that lncRNAs may function as competing endogenous RNAs (ceRNAs) to sponge microRNAs (miRNAs), thereby regulating miRNA targets (17).

Hedgehog signaling in mature osteoblasts is responsible for the pathogenesis of osteoblastic osteosarcoma through the overexpression of yes-associated protein 1 and lncRNA H19 (9). The lncRNA small nucleolar RNA host gene 1 (SNHG12) could upregulate the expression of Notch2 to promote osteosarcoma metastasis by acting as a sponge of miR-195-5p (10). The lncRNA focally amplified lncRNA on chromosome 1 has been shown to be a negative prognostic biomarker and exhibit an important pro-oncogenic effect on osteosarcoma progression by targeting epithelial-mesenchymal transition (11). Decreased expression levels of the lncRNA maternally expressed 3 (MEG3) has been reported to be an independent predictor of shorter overall survival in patients with osteosarcoma, suggesting that MEG3 lncRNA may be a useful biomarker for the prognosis of osteosarcoma (12). In addition, lncRNA MEG3 could promote osteosarcoma cell growth and metastasis in vitro by acting as a sponge of miR-127 (13). In a two-stage and case-control study of 900 osteosarcoma cases and 900 controls among the Chinese population, those with the rs7958904 CC genotype had a significantly lower expression level of the lncRNA HOX transcript antisense RNA than those with other genotypes, as well as a lower risk of osteosarcoma (14). The upregulation of lncRNA taurine upregulated 1 was also revealed to be strongly associated with poor prognosis and was an independent prognostic indicator for overall progression-free survival in osteosarcoma (15). Patients with osteosarcoma and high expression levels of lncRNA HOXA distal transcript antisense RNA had poorer overall survival than those with low expression levels (16). The majority of prognostic miRNAs were located at 14q32 in a study by Kelly et al (18); their gene targets displayed deregulation patterns associated with outcomes based on their developed miRNA models, with independent predictive values for recurrence and overall survival from formalin-fixed, paraffin-embedded human osteosarcoma diagnostic biopsy specimens. However, the roles of the lncRNAs were not investigated. In the present study, the GSE3905 dataset was downloaded to identify prognosis-associated RNAs for ceRNA regulatory network construction. In addition, the important pathways were determined using Gene Set Enrichment Analysis (GSEA). Kaplan-Meier survival analysis was conducted to determine the associations between the expression levels and survival prognosis. Furthermore, another independent miRNA profile, GSE79181, was downloaded from Gene Expression Omnibus (GEO) for validation (Fig. 1).

 Materials and methods Microarray data and preprocessing

The expression profile of GSE39058 (18), containing two subdatasets (GSE39040 and GSE39055), was downloaded from GEO (www.ncbi.nlm.nih.gov/geo/). The small non-coding RNA (miRNA) profile of GSE39040, containing 65 human osteosarcoma biopsy specimens, was generated using the platform GPL15762 (Illumina Human v2 MicroRNA Expression BeadChip; Illumina, Inc., San Diego, CA, USA) and the expression profile of GSE39055 containing 37 human osteosarcoma biopsy specimens was based on the platform of GPL14951 (Illumina HumanHT-12 WG-DASL V4.0 R2 Expression BeadChip); Illumina, Inc.). The paired miRNA and mRNA profiles of 37 biopsy specimens were used for further analysis. The probe names in the gene expression profiles were converted into gene names. In addition, the mean value of a gene symbol matched with multiple probes was calculated using the function in R software (version 3.4.1) (19). Subsequently, the microarray data were normalized using the preprocessCore package in R software (version 3.4.1) (19). The Linear Models for Microarray Analysis (Limma) package (version 3.32.5) in R software (20) was used for log2 conversion and the sequencing data were normalized using the quantile method (21).

 Differential expression analysis and hierarchical clustering

The mRNAs and lncRNAs in the GSE39055 dataset were annotated using the Transcript ID, RefSeq ID and location, supported by the GPL14951 (Illumina HumanHT-12 WG-DASL V4.0; Illumina, Inc.) platform and HUGO Gene Nomenclature Committee (www.genenames.org), which has 3,859 lncRNAs and 19,180 protein-coding genes (22). The 37 biopsy specimens of the paired miRNA and mRNA profile were divided into recurrent (n=18) and non-recurrent (n=19) osteosarcoma samples based on the clinical information. The Limma package (20) was used to perform differential expression analysis between recurrent and non-recurrent osteosarcoma samples. The P-values were adjusted for false discovery rate (FDR) using the Benjamini-Hochberg method (23). The thresholds for screening differentially expressed mRNAs, miRNAs and lncRNA were set as FDR<0.05 and |log2 fold-change (FC)|>0.585. In addition, the expression values of screened differentially expressed RNAs were hierarchically clustered using the pheatmap package (version 1.0.8) (24) in R based on the encyclopedia of distances (25) to observe the differences in expression levels spontaneously.

 Screening prognosis-associated RNAs

The expression levels of differentially expressed RNAs and the prognostic information of each sample were collected. The survival package (26) in R was used to identify the prognosis-associated mRNAs, miRNAs and lncRNAs by Cox regression univariate analysis, and P<0.05 was used as the cut-off criterion. Subsequently, the Database for Annotation, Visualization and Integrated Discovery (version 6.8) (27) was used to perform Gene Ontology (GO) functional [biology process (BP), molecular function (MF) and cellular component (CC)] and Kyoto Encyclopedia of Genes and Genomes (KEGG) pathway enrichment analysis for prognosis-associated mRNAs, with a threshold of P<0.05.

 ceRNA regulatory network construction and analysis

The identified prognosis-associated lncRNAs and miRNAs were used to select lncRNA-miRNA regulatory relationships from the associations in miRcode (version 11; www.mircode.org/) and starBase (version 2.0; starbase.sysu.edu.cn/index.php) (28). In addition, the miRNA-mRNA regulatory associations were selected based on the regulatory information in TargetScan (release 7.1; www.targetscan.org/vert_71/) (29), miRTarBase (release 6.0; mirtarbase.mbc.nctu.edu.tw) (30), miRanda (version 2.0; www.microrna.org/microrna/home.do) (31) and miRBase (release 21; www.mirbase.org/) (32). The ceRNA regulatory network was constructed by integrating the lncRNA-miRNA and miRNA-mRNA regulatory relationships via Cytoscape (version 3.4; www.cytoscape.org/) (33).

The genes in the constructed ceRNA regulatory network were used for GSEA (2017 update; software.broadinstitute.org/gsea/index.jsp) to identify the important KEGG pathways (34). Subsequently, the KEGG pathways directly associated with osteosarcoma were screened by searching the term ‘osteosarcoma’ in the Comparative Toxicogenomics Database (CTD; 2017 update; ctd.mdibl.org/) (35).

 Validation of prognosis-associated RNAs

The important RNAs participating in the KEGG pathways directly associated with osteosarcoma were used to classify the samples into high- and low-expression groups based on the expression levels. Kaplan-Meier survival analysis, conducted with R software (version 3.4.1), was used to determine the associations between the expression levels and survival prognosis. Another independent miRNA profile, GSE79181, containing 25 osteosarcoma samples was downloaded from GEO (www.ncbi.nlm.nih.gov/geo/) to validate the reliability of the important molecules.

 Results DEGs screening and hierarchical clustering

Following data preprocessing, 15 differentially expressed lncRNAs (8 down-regulated and 7 upregulated lncRNAs) and 837 differentially expressed mRNAs (362 downregulated and 475 upregulated mRNAs) were identified in the GSE30955 profile. In addition, 32 differentially expressed miRNAs, including 13 down-regulated and 19 upregulated miRNAs, were revealed in the GSE39040 profile. The expression values of differentially expressed miRNAs, lncRNAs, and mRNAs were hierarchically clustered, and the color contrast indicated a significant difference in the expression levels between the recurrent and non-recurrent osteosarcoma samples (Fig. 2).

 Prognosis-associated mRNAs, miRNAs and lncRNAs

Among the selected differentially expressed RNAs, 417 RNAs (5 lncRNAs, 19 miRNAs and 393 mRNAs) were revealed to be associated with prognosis. The functional enrichment analysis for prognosis-associated mRNAs demonstrated that these 393 mRNAs were markedly involved with 10, 8 and 2 GO terms in the BP, CC and MF categories, respectively (Table I). These functions were mainly associated with ‘chromatin organization’ and ‘transcription regulator’. Furthermore, 9 KEGG pathways were enriched for prognosis-associated mRNAs, such as ‘transforming growth factor-β signaling pathway’ (P= 0.0258), ‘Wnt signaling pathway’ (P= 0.0366) and ‘Oxidative phosphorylation’ (P=0.0493; Table II).

 ceRNA regulatory network construction and analysis

A total of 65 lncRNA-miRNA regulatory relationships associated with the 5 prognosis-associated lncRNAs were identified in miRcode and starBase. In addition, 7 lncRNA-miRNA regulatory relationships associated with the 19 prognosis-associated miRNAs were further screened, which involved three lncRNAs [long intergenic non-protein coding RNA 28 (LINC00028), LINC00323 and SNHG1] and four miRNAs (hsa-miR-124, hsa-miR-139-5p, hsa-miR-7, and hsa-miR-206). Furthermore, 172 miRNA-mRNA regulatory relationships were obtained for these four miRNAs using TargetScan, miRBase, miRanda, and miRTarBase. Finally, the 7 lncRNA-miRNA and 172 miRNA-mRNA regulatory relationships were used to construct the ceRNA regulatory network, which included 3 lncRNAs, 4 miRNAs and 143 mRNAs (Fig. 3).

The GSEA for the ceRNA regulatory network revealed that three KEGG pathways [‘Mitogen-activated protein kinase (MAPK) signaling pathway’, ‘Chemokine signaling pathway’ and ‘Spliceosome’] were markedly associated with osteosarcoma (Table III). Furthermore, ‘MAPK signaling pathway’ involving three genes [Ras-related protein Rap-1b (RAP1B), activating transcription factor 2 (ATF2) and protein phosphatase Mg2+/Mn2+ dependent 1B (PPM1B)] was also demonstrated to be associated with osteosarcoma by searching the CTD. Based on the ceRNA regulatory network, two miRNAs (hsa-miR-124 and hsa-miR-7) and three lncRNAs (LINC00028, LINC00323, and SNHG1) could regulate the three mRNAs (RAP1B, ATF2, and PPM1B) that participate in the MAPK signaling pathway (Fig. 4).

 Validation of prognosis-associated RNAs

The samples were classified into high- and low-expression groups based on the expression levels of three mRNAs (RAP1B, ATF2, and PPM1B), two miRNAs (hsa-miR-124 and hsa-miR-7), and three lncRNAs (LINC00028, LINC00323, and SNHG1). Kaplan-Meier survival analysis revealed that samples with lower expression levels of LINC00323 and SNHG1 had better prognosis, and increased expression levels of LINC00028, hsa-miR-124, and hsa-miR-7 were associated with better prognosis. The overexpression of RAP1B, ATF2 and PPM1B was also positively associated with osteosarcoma recurrence (Fig. 5).

The expression level of hsa-miR-7 was significantly down-regulated in the osteosarcoma recurrent samples (P=0.02346) of another independent miRNA profile, GSE79181 (Fig. 6A). Kaplan-Meier survival analysis revealed that hsa-miR-124 (P=0.04681) and hsa-miR-7 (P=0.04853) were markedly associated with osteosarcoma recurrence (Fig. 6) using GSE79181. These results were in accordance with the GSE39040 analysis results.

 Discussion

Understanding the mechanisms underlying osteosarcoma pathogenesis is necessary to determine novel biomarkers for the development of novel treatment strategies, particularly for metastases, as the prognosis of osteosarcoma is still poor (36). A number of previous studies have suggested that the dysregulated expression of lncRNAs could independently predict the outcomes of patients with osteosarcoma (37,38). In the present study, 417 RNAs (5 lncRNAs, 19 miRNAs and 393 mRNAs) were revealed to be associated with prognosis among the total number of differentially expressed RNAs identified. Subsequently, the prognosis-associated RNAs were used to construct a ceRNA regulatory network for GSEA. The GSEA for the ceRNA regulatory network revealed that ‘MAPK signaling pathway’, ‘Chemokine signaling pathway’ and ‘Spliceosome’ were markedly associated with osteosarcoma. In addition, the MAPK signaling pathway was also demonstrated to be associated with osteosarcoma via the CTD.

MAPK is an insulin-mitogen activated protein (Ser/Thr) kinase, and activated MAPK can transmit extracellular signals to regulate cell growth, migration and apoptosis (39). Many studies have reported the significant roles of the MAPK signaling pathway in osteosarcoma (40,41). The expression level of S100 calcium binding protein A9 (S100A9) was increased in human osteosarcoma tissues and associated with survival rate, and the downregulation of S100A9 inhibited osteosarcoma cell proliferation through the inactivation of the MAPK and NF-κB signaling pathways (42). Increased expression levels of phospholipase A2 group XVI could enhance osteosarcoma metastasis by activating the MAPK signaling pathway (43). Delphinidin treatment could inhibit cell migration and prevented epithelial-to-mesenchymal transition via the MAPK signaling pathway in osteosarcoma cell lines (44). Onzin overexpression could also contribute to more aggressive osteosarcoma metastatic phenotypes through the C-X-C motif chemokine 5-MAPK signaling pathway (45).

In the constructed ceRNA regulatory network, three lncRNAs (LINC00028, LINC00323 and SNHG1) and two miRNAs (hsa-miR-124 and hsa-miR-7) regulating three mRNAs (RAP1B, ATF2 and PPM1B) participated in the MAPK signaling pathway. In addition, samples with increased expression levels of hsa-miR-124 and hsa-miR-7 had better prognosis. Expression levels of miR-124 in the metastatic osteosarcoma tissues were revealed to be significantly lower than those observed in non-metastatic tissues, and miR-124 could function as a tumor suppressor miRNA in osteosarcoma through the inhibition of Rac family small guanosine triphosphatase 1 expression (46). In addition, miR-7 expression was decreased in osteosar-coma tissues and low miR-7 expression was associated with poor prognosis (47). Furthermore, the roles of hsa-miR-124 and hsa-miR-7 in osteosarcoma recurrence were also validated using another independent miRNA profile in the present study.

Kaplan-Meier survival analysis revealed that samples with lower expression levels of LINC00323 and SNHG1 had better prognosis, which were consistent with a study by Wang et al (48), which demonstrated that SNHG1 was upregulated in osteosarcoma tissues and cell lines, and high SNHG1 expression predicted poor overall survival in patients with osteosarcoma; these findings were consistent with the present results. Wang et al (48) reported that SNHG1 increased human nin one binding protein (an oncogene) by acting as a ceRNA to sponge miR-326, promoting cell growth, migration and invasion in osteosarcoma, while the present study revealed that SNHG1 increased ATF2, PPM1B and RAP1B by sponging miR-124 and miR-7 as a ceRNA. In addition, high SNHG1 expression was also positively associated with the metastasis of osteosarcoma, and miR-577 could act as a ceRNA of SNHG1 by targeting WNT2B which served an oncogenic role in osteosarcoma cells by activating the Wnt/β-catenin signaling pathway (49). Therefore, the present study showed that the SNHG1/miR-124 and miR-7/RAP1B, ATF2 and PPM1B/MAPK signaling pathway regulatory network may provide a potential novel therapeutic strategy for osteosarcoma treatment.

In conclusion, three lncRNAs (LINC00028, LINC00323 and SNHG1) and two miRNAs (hsa-miR-124 and hsa-miR-7) regulating three mRNAs (RAP1B, ATF2 and PPM1B) participated in MAPK signaling pathway, as identified by the ceRNA regulatory network, and were associated with osteosarcoma recurrence. However, the prognosis prediction value based on the expressions of these important RNAs still requires further confirmation based on clinical experiments in an independent cohort of patients with non-recurrent and recurrent osteosarcoma.

 Funding

No funding was received.

 Availability of data and materials

The datasets used and/or analyzed during the present study are available from the corresponding author on reasonable request.

 Authors' contributions

SZ and HF contributed to the study design. SZ, LD, XL and HF searched for and downloaded the gene expression profiles from the Gene Expression Omnibus database. SZ, LD, XL and HF performed the analysis and interpretation of the microarray dataset. SZ and HF critically revised the manuscript for important intellectual content. All authors read and approved the final manuscript.

 Ethics approval and consent to participate

Not applicable.

 Patient consent for publication

Not applicable.

 Competing interests

The authors declare that they have no competing interests.

 Acknowledgments

Not applicable.

 References GellerDSGorlickROsteosarcoma: A review of diagnosis, management, and treatment strategiesClin Adv Hematol Oncol87057182010 HeHNiJHuangJMolecular mechanisms of chemoresistance in osteosarcoma (Review)Oncol Lett713521362201410.3892/ol.2014.1935247651373997672 LiuJJLiuSWangJGZhuWHuaYQSunWCaiZDTelangiectatic osteosarcoma: A review of literatureOnco Targets Ther65936022013237450513671797 McQueenPGhaffarSGuoYRubinEMZiXHoangBHThe Wnt signaling pathway: Implications for therapy in osteosarcomaExpert Rev Anticancer Ther1112231232201110.1586/era.11.9421916576 LinCHJiTChenCFHoangBHWnt signaling in osteosarcomaAdv Exp Med Biol8043345201410.1007/978-3-319-04843-7_224924167 McManusMMWeissKRHughesDPUnderstanding the role of Notch in osteosarcomaAdv Exp Med Biol8046792201410.1007/978-3-319-04843-7_424924169 GrazianoACCardileVAvolaRVicarioNParentiCSalvatorelliLMagroGParentiRWilms' tumor gene 1 silencing inhibits proliferation of human osteosarcoma MG-63 cell line by cell cycle arrest and apoptosis activationOncotarget81391713931201710.18632/oncotarget.14715281071965355150 ZhangFYanTGuoWSunKWangSBaoXLiuKZhengBZhangHRenTNovel oncogene COPS3 interacts with beclin1 and Raf-1 to regulate metastasis of osteosarcoma through autophagyJ Exp Clin Cancer Res37135201810.1186/s13046-018-0791-6299701156029018 ChanLHWangWYeungWDengYYuanPMakKKHedgehog signaling induces osteosarcoma development through Yap1 and H19 overexpressionOncogene3348574866201410.1038/onc.2013.433 ZhouSYuLXiongMDaiGLncRNA SNHG12 promotes tumorigenesis and metastasis in osteosarcoma by upregulating Notch2 by sponging miR-195-5pBiochem Biophys Res Commun49518221832201810.1016/j.bbrc.2017.12.047 WangYZhaoZZhangSLiZLiDYangSZhangHZengXLiuJLncRNA FAL1 is a negative prognostic biomarker and exhibits pro-oncogenic function in osteosarcomaJ Cell Biochem11984818489201810.1002/jcb.2707429987852 TianZZGuoXJZhaoYMFangYDecreased expression of long non-coding RNA MEG3 acts as a potential predictor biomarker in progression and poor prognosis of osteosarcomaInt J Clin Exp Pathol815138151422015 WangYKongDKnockdown of lncRNA MEG3 inhibits viability, migration, and invasion and promotes apoptosis by sponging miR-127 in osteosarcoma cellJ Cell Biochem119669679201810.1002/jcb.26230 ZhouQChenFFeiZZhaoJLiangYPanWLiuXZhengDGenetic variants of lncRNA HOTAIR contribute to the risk of osteosarcomaOncotarget719928199342016269673894991428 MaBLiMZhangLHuangMLeiJBFuGHLiuCXLaiQWChenQQWangYLUpregulation of long non-coding RNA TUG1 correlates with poor prognosis and disease status in osteosarcomaTumour Biol3744454455201610.1007/s13277-015-4301-6 LiFCaoLHangDWangFWangQLong non-coding RNA HOTTIP is up-regulated and associated with poor prognosis in patients with osteosarcomaInt J Clin Exp Pathol811414114202015266178684637684 SalmenaLPolisenoLTayYKatsLPandolfiPPA ceRNA hypothesis: The rosetta stone of a hidden RNA languageCell146353358201110.1016/j.cell.2011.07.014218021303235919 KellyADHaibe-KainsBJanewayKAHillKEHoweEGoldsmithJKurekKPerez-AtaydeARFrancoeurNFanJBMicroRNA paraffin-based studies in osteosarcoma reveal reproducible independent prognostic profiles at 14q32Genome Med52201310.1186/gm406233394623706900 BolstadBPreprocessCore: A collection of pre-processing functionsR package version 12013 SmythGKLimma: Linear models for microarray data Bioinformatics computational biology Solutions Using R BioconductorSpringer3974202005 RaoYLeeYJarjouraDRuppertASLiuCGHsuJCHaganJPA comparison of normalization techniques for microRNA microarray dataStat Appl Genet Mol Biol7Article22200810.2202/1544-6115.128718673291 YatesBBraschiBGrayKASealRLTweedieSBrufordEAGenenames.org: The HGNC and VGNC resources in 2017Nucleic Acids Res45D619D625201710.1093/nar/gkw1033 BenjaminiYHochbergYControlling the false discovery rate: A practical and powerful approach to multiple testingJ Royal Stat Soc Series B572893001995 WangLCaoCMaQZengQWangHChengZZhuGQiJMaHNianHWangYRNA-seq analyses of multiple meristems of soybean: Novel and alternative transcripts, evolutionary and functional implicationsBMC Plant Biol14169201410.1186/1471-2229-14-169249395564070088 GoslingCEncyclopedia of distancesReference Rev2415832009 WangPWangYHangBZouXMaoJHA novel gene expression-based prognostic scoring system to predict survival in gastric cancerOncotarget755343553512016274193735342421 Huang DaWShermanBTLempickiRASystematic and integrative analysis of large gene lists using DAVID bioinformatics resourcesNature Protoc44457200910.1038/nprot.2008.211 LiJHLiuSZhouHQuLHYangJHstarBase v2.0: Decoding miRNA-ceRNA, miRNA-ncRNA and protein-RNA interaction networks from large-scale CLIP-Seq dataNucleic Acids Res42D92D97201410.1093/nar/gkt1248 AgarwalVBellGWNamJWBartelDPPredicting effective microRNA target sites in mammalian mRNAsElife1242015 LeeYYKimTJKimJYChoiCHDoIGSongSYSohnIJungSHBaeDSLeeJWKimBGGenetic profiling to predict recurrence of early cervical cancerGynecol Oncol131650654201310.1016/j.ygyno.2013.10.00324145113 JeggariAMarksDSLarssonEmiRcode: A map of putative microRNA target sites in the long non-coding transcriptomeBioinformatics2820622063201210.1093/bioinformatics/bts344227187873400968 KozomaraAGriffiths-JonesSmiRBase: Annotating high confidence microRNAs using deep sequencing dataNucleic Acids Res42D68D73201410.1093/nar/gkt11813965103 ShannonPMarkielAOzierOBaligaNSWangJTRamageDAminNSchwikowskiBIdekerTCytoscape: A software environment for integrated models of biomolecular interaction networksGenome Res1324982504200310.1101/gr.123930314597658403769 SubramanianATamayoPMoothaVKMukherjeeSEbertBLGilletteMAPaulovichAPomeroySLGolubTRLanderESMesirovJPGene set enrichment analysis: A knowledge-based approach for interpreting genome-wide expression profilesProc Natl Acad Sci USA1021554515550200510.1073/pnas.0506580102161995171239896 DavisAPGrondinCJJohnsonRJSciakyDKingBLMcMorranRWiegersJWiegersTCMattinglyCJThe comparative toxicogenomics database: Update 2017Nucleic Acids Res45D972D978201710.1093/nar/gkw8385210612 HeJPHaoYWangXLYangXJShaoJFGuoFJFengJXReview of the molecular pathogenesis of osteosarcomaAsian Pac J Cancer Prev1559675976201410.7314/APJCP.2014.15.15.596725124559 WangBSuYYangQLvDZhangWTangKWangHZhangRLiuYOverexpression of long non-coding RNA HOTAIR promotes tumor growth and metastasis in human osteosarcomaMol Cells38432440201510.14348/molcells.2015.2327257287534443285 ZhangQGengPLYinPWangXLJiaJPYaoJDown-regulation of long non-coding RNA TUG1 inhibits osteo-sarcoma cell proliferation and promotes apoptosisAsian Pac J Cancer Prev1423112325201310.7314/APJCP.2013.14.4.2311 SuiXKongNYeLHanWZhouJZhangQHeCPanHp38 and JNK MAPK pathways control the balance of apoptosis and autophagy in response to chemotherapeutic agentsCancer Lett344174179201410.1016/j.canlet.2013.11.019 ChenHJLinCMLeeCYShihNCPengSFTsuzukiMAmagayaSHuangWWYangJSKaempferol suppresses cell metastasis via inhibition of the ERK-p38-JNK and AP-1 signaling pathways in U-2 OS human osteosarcoma cellsOncol Rep30925932201310.3892/or.2013.249023708932 ChengDDZhuBLiSJYuanTYangQCFanCYDown-regulation of RPS9 inhibits osteosarcoma cell growth through inactivation of MAPK signaling pathwayJ Cancer827202728201710.7150/jca.19130289288615604204 ChengSZhangXHuangNQiuQJinYJiangDDown-regulation of S100A9 inhibits osteosarcoma cell growth through inactivating MAPK and NF-κB signaling pathwaysBMC Cancer16253201610.1186/s12885-016-2294-1 LiLLiangSWasylishenARZhangYYangXZhouBShanLHanXMuTWangGXiongSPLA2G16 promotes osteosarcoma metastasis and drug resistance via the MAPK pathwayOncotarget718021180352016269338044951268 KangHMParkBSKangHKParkHRYuSBKimIRDelphinidin induces apoptosis and inhibits epithelial-to-mesenchymal transition via the ERK/p38 MAPK-signaling pathway in human osteosarcoma cell linesEnviron Toxicol33640649201810.1002/tox.22548294513515969316 ZhangYHuQLiGLiLLiangSZhangYLiuJFanZLiLZhouBONZIN upregulation by mutant p53 contributes to osteosarcoma metastasis through the CXCL5-MAPK signaling pathwayCell Physiol Biochem4810991111201810.1159/00049197630041188 GengSZhangXChenJLiuXZhangHXuXMaYLiBZhangYBiZYangCThe tumor suppressor role of miR-124 in osteosarcomaPLoS one9e91566201410.1371/journal.pone.0091566249719024074038 ShouyingLChangxiZChangbaoZQiuheSMingWYeLHuaijieALow-expression of miR-7 promotes cell proliferation and exhibits prognostic value in osteosarcoma patientsInt J Clin Exp Pathol10903590412017 WangJCaoLWuJWangQLong non-coding RNA SNHG1 regulates NOB1 expression by sponging miR-326 and promotes tumorigenesis in osteosarcomaInt J Oncol5277882018 JiangZJiangCFangJUp-regulated lnc-SNHG1 contributes to osteosarcoma progression through sequestration of miR-577 and activation of WNT2B/Wnt/β-catenin pathwayBiochem Biophys Res Commun495238245201810.1016/j.bbrc.2017.11.012

Flow chart of the analysis performed in the present study. GSE39058 was downloaded in order to identify prognosis-associated lncRNAs, miRNAs and mRNAs for ceRNA regulatory network construction. Important pathways were searched using GSEA. Kaplan-Meier survival analysis was conducted to determine the associations between the expression levels and survival prognosis. In addition, another independent dataset GSE79181 was also downloaded for validation. GSEA, Gene Set Enrichment Analysis; lncRNA, long non-coding RNA; miRNA, microRNA; ceRNA, competing endogenous RNA; CTD, Comparative Toxicogenomics Database.

Hierarchical clustering analysis of screened differentially expressed RNAs. The color contrast indicates the marked differences in the expression levels between the recurrent and non-recurrent osteosarcoma samples. (A) Heatmap for GSE39040 using the differentially expressed miRNAs. (B) Heatmap for GSE39055 using the differentially expressed lncRNAs and mRNAs separately. lncRNA, long non-coding RNA; miRNA, microRNA.

Constructed ceRNA regulatory network using prognosis-associated RNAs. The ceRNA regulatory network included 7 lncRNA-miRNA and 172 miRNA-mRNA regulatory associations, which involved 3 lncRNAs, 4 miRNAs and 143 mRNAs. Red nodes indicate the upregulated RNAs, and green indicate the downregulated RNAs. lncRNA, long non-coding RNA; miRNA/miR, microRNA; ceRNA, competing endogenous RNA; FC, fold change.

ceRNA regulatory network associated with the MAPK signaling pathway. Two miRNAs (hsa-miR-124 and hsa-miR-7) and three lncRNAs (LINC00028, LINC00323, and SNHG1) could regulate three mRNAs (RAP1B, ATF2 and PPM1B) that participate in the MAPK signaling pathway. lncRNA, long non-coding RNA; miRNA/miR, microRNA; ceRNA, competing endogenous RNA; LINC, long intergenic non-protein coding RNA; SNHG1, small nucleolar RNA host gene 1; RAP1B, Ras-related protein Rap-1b; ATF2, activating transcription factor 2; PPM1B, protein phosphatase Mg2+/Mn2+ dependent 1B; MAPK, mitogen-activated protein kinase.

Kaplan-Meier survival analysis for lncRNAs, miRNAs and mRNAs associated with the mitogen-activated protein kinase signaling pathway. (A) Lower expression levels of LINC00323 and SNHG1, and increased expression levels of LINC00028 were associated with better prognosis. (B) Samples with increased expression levels of hsa-miR-124 and hsa-miR-7 had better prognosis. (C) The overexpression of RAP1B, ATF2, and PPM1B was positively associated with osteosarcoma recurrence. lncRNA, long non-coding RNA; miRNA/miR, microRNA; ceRNA, competing endogenous RNA; LINC, long intergenic non-protein coding RNA; SNHG1, small nucleolar RNA host gene 1; RAP1B, Ras-related protein Rap-1b; ATF2, activating transcription factor 2; PPM1B, protein phosphatase Mg2+/Mn2+ dependent 1B.

Kaplan-Meier survival analysis of hsa-miR-124 and hsa-miR-7 with osteosarcoma recurrence using GSE79181 for validation. (A) hsa-miR-7 was significantly downregulated in the osteosarcoma recurrent samples. (B) hsa-miR-124 was markedly associated with osteosarcoma recurrence. Data are presented as the mean ± standard deviation. *P<0.05 vs. non-recurrence. miR, microRNA.

Gene Ontology functional enrichment analysis in terms of biology process, cellular component and molecular function for the prognosis-associated mRNAs.

Category Term Count P-Value Genes
BP GO:0043009~chordate embryonic development 16 0.0033 MAFG, NBN, SYVN1, ADA, HOXA2, HOXA3, HOXC9, MEOX2, HOXA5, HOXB8, NCOA6, APBA3, DAD1, MAPK8IP3, ROR2, MED1
GO:0048568~embryonic organ development 10 0.0086 HOXA2, HOXC9, HOXA3, HOXB8, HOXA5, MYO7A, ROR2, VAX2, ADA, MED1
GO:0003002~regionalization 10 0.0195 HOXA2, HOXC9, HOXA3, BTG2, MEOX2, HOXB8, HOXA5, ROR2, VAX2, LRP5
GO:0051276~chromosome organization 18 0.0208 MEAF6, NBN, HIST1H1D, HIST1H2BF, HP1BP3, HIST1H1A, EZH2, CENPE, TSPYL1, CHMP1A, EYA2, HDAC2, NCAPG, PRDM5, CHD1,TLK1, RNF40, KDM5D
GO:0001501~skeletal system development 13 0.0296 ESRRA, AEBP1, BMP1, TUFT1, ATP6V1B1, HOXA2, HOXA3, HOXC9, HOXB8, HOXA5, ETS2, COMP, ROR2
GO:0006351~transcription, DNA-dependent 12 0.0367 TAF11, TAF1B, MED7, CCNK, TAF1A, ETS1, ANG, NCOA6, NFIX, POLR2D, MED1, TAF1L
GO:0032774~RNA biosynthetic process 12 0.0399 TAF11, TAF1B, MED7, CCNK, TAF1A, ETS1, ANG, NCOA6, NFIX, POLR2D, MED1, TAF1L
GO:0045597~positive regulation of cell differentiation 10 0.0448 IKZF1, ETS1, VHL, MAPT, KITLG, SEMA4D, SOCS5, NTN1, BOC, ADA
GO:0006325~chromatin organization 14 0.0458 MEAF6, HIST1H1D, HP1BP3, HIST1H2BF, HIST1H1A, EZH2, TSPYL1, EYA2, HDAC2, PRDM5, CHD1, TLK1, RNF40, KDM5D
GO:0048598~embryonic morphogenesis 12 0.0495 HOXA2, EYA2, HOXC9, HOXA3, HOXB8, HOXA5, MYO7A, CRABP2, ROR2, VAX2, MED1, LRP5
CC GO:0044463~cell projection part 15 0.0002 PLD2, BBS5, PARD3, TTLL6, KLC2, GAS8, ADA, SPAG16, PCSK1, CC2D2A, CYBRD1, MAPK8IP3, PEBP1, SLC5A6, CNTNAP1
GO:0042995~cell projection 26 0.0017 BBS5, PARD3, MYO7A, CLSTN1, TTLL6, KLC2, ADA, SPAG16, PCSK1, ANG, MAPT, SLC4A7, CNTNAP1, PLD2, DBNL, ARHGEF7, NPY1R, GAS8, THY1, GRM2, CC2D2A, CYBRD1, PEBP1, MAPK8IP3, SLC5A6, GAP43
GO:0005694~chromosome 19 0.0030 TCP1, NBN, HIST1H1D, PRPF4B, IKZF1, HIST1H2BF, PPP2R5A, HP1BP3, HIST1H1A, AKAP8L, CENPE, TRIM66, CHMP1A, HDAC2, RIF1, NCAPG, PPP2CB, CHD1, RNF40
GO:0044427~chromosomal part 15 0.0158 TCP1, NBN, HIST1H1D, IKZF1, HIST1H2BF, HP1BP3, PPP2R5A, HIST1H1A, CENPE, TRIM66,HDAC2, RIF1, NCAPG, PPP2CB, CHD1
GO:0005578~proteinaceous extracellular matrix 13 0.0195 ASPN, COL4A2, HMCN1, ADAMTSL2, LAMC3, ANG, FBLN2, COMP, TIMP4, MMP3, NTN1, MFAP5, SPON1
GO:0043005~neuron projection 13 0.0306 ARHGEF7, NPY1R, KLC2, ADA, THY1, PCSK1, GRM2, ANG, MAPT, MAPK8IP3, PEBP1, CNTNAP1, GAP43
GO:0031012~extracellular matrix 13 0.0322 ASPN, COL4A2, HMCN1, ADAMTSL2, LAMC3, ANG, FBLN2, COMP, TIMP4, MMP3, NTN1, MFAP5, SPON1
GO:0005730~nucleolus 21 0.0440 UTP23, MEAF6, PNMA3, NBN, SYVN1, TAF1A, TSR1, HNRNPA1, PNN, TSPYL1, NOP14, STAT6, DIS3, KRT17, ANG, RPL9, NCOA6, PQBP1, LIAS, HNRNPH1, ARL4A
MF GO:0043565~sequence-specific DNA binding 22 0.0116 MAFG, TSHZ3, ERF, ESRRA, IKZF1, TFE3, VAX2, NFIX, ATF2, STAT6, ATF5, HOXA2, HDAC2, HOXA3, HOXC9, MEOX2, HOXA5, ETS1, HOXB8, ETS2, PRDM5, CREB3L4
GO:0030528~transcription regulator activity 41 0.0440 TAF1B, TSHZ3, AEBP1, TAF1A, TFE3, E2F8, EZH2, CNOT3, NFIX, MXI1, ATF2, ZNF75D, KCNIP3, STAT6, TRIM66, HOXA2, HOXC9, HOXA3, HOXA5, PQBP1, CREB3L4, TAF1L, MAFG, ERF, ESRRA, IKZF1, VAX2, AFF1, MSRB2, TAF11, MED7, ATF5, HDAC2, MEOX2, ETS1, HOXB8, ETS2, PRDM5, NCOA6, ID4, MED1

GO, Gene Ontology; BP, biology process; CC, cellular component; MF, molecular function.

Nine enriched Kyoto Encyclopedia of Genes and Genomes pathways for prognosis-associated mRNAs.

KEGG Pathway Count P-Value Genes
hsa00520:Amino sugar and nucleotide sugar metabolism 3 0.0224 GMPPB, PGM3, FUK
hsa00603:Glycosphingolipid biosynthesis 2 0.0249 GBGT1, B3GALNT1
hsa04350:TGF-β signaling pathway 4 0.0258 AMHR2, COMP, PPP2CB, ID4
hsa04144:Endocytosis 6 0.0316 PRKCZ, PLD2, PARD3, RAB11FIP3, ACAP2, AGAP1
hsa04310:Wnt signaling pathway 5 0.0366 PPP2R5A, PPP2CB, PRICKLE2, PLCB1, LRP5
hsa03020:RNA polymerase 2 0.0437 POLR2D, POLR2J2
hsa04070:Phosphatidylinositol signaling system 3 0.0445 DGKQ, PLCD4, PLCB1
hsa03040:Spliceosome 4 0.0473 PQBP1, MAGOHB, THOC2, HNRNPA1
hsa00190:Oxidative phosphorylation 4 0.0493 NDUFB4, ATP6V1E2, ATP6V1B1, COX17

KEGG, Kyoto Encyclopedia of Genes and Genomes; TGF-β, transforming growth factor-β.

Gene Set Enrichment Analysis for the competing endogenous RNA regulatory network.

Name Enrichment score Normalized enrichment score Nominal P-value Genes
MAPK_SIGNALING_PATHWAY 0.5145 1.1238 0.0324 RAP1B, ATF2, PPM1B
CHEMOKINE_SIGNALING_PATHWAY 0.4710 1.0218 0.0423 RAP1B, PLCB1, PRKCZ
SPLICEOSOME 0.3188 0.6929 0.0485 MAGOHB, PQBP1

MAPK, mitogen-activated protein kinase; RAP1B, Ras-related protein Rap-1b; ATF2, activating transcription factor 2; PPM1B, protein phosphatase Mg2+/Mn2+ dependent 1B; PLCB1, phospholipase Cβ1; PRKCZ, protein kinase Cζ; MAGOHB, Mago homolog B exon junction complex subunit; PQBP1, polyglutamine binding protein 1.