miRNA and mRNA expression profiles for osteosarcoma prior to and following treatment were obtained from the GEO database (accession nos. GSE39040 and GSE39055, respectively; https://www.ncbi.nlm.nih.gov/geo/) (15). Probes were mapped to genes; a gene with numerous probes was represented by the mean of the expression values of the probes. Matching miRNA and mRNA expression profiles were obtained for further analysis. Clinical data for the patients with osteosarcoma were also obtained for the two datasets, including age, sex, recurrence status and other clinical information.

Human miRNA-mRNA target data was collected from TarBase (16), miRTarBase (17) and miRecords (18), which manually curate data regarding high-quality, experimentally validated miRNA-mRNA target interactions from published experiments. By integrating data from these three databases, a total of 44,181 non-redundant miRNA-mRNA target interactions were obtained.

The OSceNet was constructed using R software (v3.2.1; http://mirrors.tuna.tsinghua.edu.cn/CRAN/) according to two principles: The two mRNAs were regulated by the same miRNA, and the two mRNAs followed the same expression pattern. Firstly, correlations between the previously validated miRNA-mRNA target pairs in the matched miRNA and mRNA expression profiles were evaluated by Pearson correlation coefficients. The pairs with a significant negative correlation were regarded as relevant miRNA-mRNA pairs in the context of osteosarcoma (R< −0.4, adjusted P<0.05). Secondly, all mRNA-mRNA pairs that were both targeted by at least one miRNA were listed as candidate ceRNA pairs. Thirdly, for each candidate ceRNA pair, the Pearson correlation coefficient between them was calculated. All candidate ceRNA pairs with R>0.4 and adjusted P<0.05 were considered ceRNA-ceRNA interactions. By assembling all the identified ceRNA pairs, the OSceNet was constructed. Within the OSceNet, each node represents an mRNA, and two nodes are connected if they are co-regulated by at least one miRNA and are co-expressed in osteosarcoma.

The networks were visualized using Cytoscape 3.3.0 (19), including the OSceNet and the modules view.

The present study assessed the impact on the recurrence-free prognosis for each module of the OSceNet using the risk score model. The risk score model was constructed by considering the power of each of the ceRNAs in the module, as evaluated by univariate Cox regression analysis, and the relative expression of each ceRNA, as follows: Risk score = ∑ i = 1 n c o e f i ∗ e x p iwhere n is the number of ceRNAs in the module, exp_i is the expression level of ceRNA i, and coef_i is the estimated regression coefficient for ceRNA i from the univariate Cox regression analysis. The median risk score value was selected as a cutoff to classify patients into high- and low-risk groups. Recurrence-free survival analyses were performed to assess the difference in recurrence probability between the high- and low-risk groups, and the statistical significance was determined by a log-rank test using the R package ‘survival’ (https://github.com/therneau/survival).

Random networks with the same architecture in OSceNet were generated using R package ‘igraph’ (http://igraph.org). Modules in OSceNet were identified using cFinder software (20). Multivariate Cox regression analyses were performed using Cox proportional hazards regression model to determine whether the prognostic models were independent of other clinical variables, adjusting for age, sex, primary disease site and chemosensitivity as covariates. The time-dependent receiver operating characteristic (ROC) curve analysis was performed using R software. All statistical analyses were performed using R software (v3.2.1), including Fisher’s exact test and Wilcoxon sum rank test. In the present study, P<0.1 was considered to indicate a statistically significant difference due to the small sample size (37 patients).

The matched miRNA and mRNA expression profiles of 37 patients with osteosarcoma were obtained from the GEO database, after which the OSceNet was constructed. Three steps were followed to generate the network: i) A total of 130 miRNA-mRNA pairs were determined between 28 miRNAs and 122 mRNAs, which were significantly negatively correlated in osteosarcoma (R< -0.4, adjusted P<0.05). ii) A total of 1,826 candidate ceRNA-ceRNA pairs were identified, which were co-regulated by at least one miRNA. iii) A total of 156 candidate ceRNA-ceRNA pairs were identified as positively correlated with each other (R>0.4, adjusted P<0.05), forming a network of 64 mRNAs (Fig. 1A).

To explore the structure and features of the OSceNet, degree distribution analysis of the nodes in the network was performed. The present study identified that most nodes had relatively few interactions, whereas a small proportion of nodes had several interactions, which fits with a power-law distribution, thus suggesting that the network was scale-free and different from randomly generated networks (Fig. 1B). Subsequently, 1,000 random networks were generated that preserved the same number of nodes and edges, and the same degree for each node as in OSceNet using the R package ‘igraph’. Comparing the clustering coefficient between the OSceNet and the random networks revealed strong, non-random competitive interactions between the OSceNet ceRNA-ceRNA pairs (Fig. 1C; P<0.001).

To analyze the competitive intensity among the ceRNA interaction pairs in OSceNet, the present study demonstrated that the nodes with a higher degree, indicating more competitors, were more strongly regulated by their interacting neighbors (R=0.57, P=0.056; Fig. 2A), which was consistent with a previous study (21). In addition, the present study investigated the effects of miRNA expression on ceRNA-ceRNA interactions. For each miRNA shared by a ceRNA pair, the tumor expression profiles were divided into high and low miRNA expression, using the median value for miRNA expression as the cutoff. By comparing the co-expression of the ceRNA pairs between these two groups of tumor samples, it was revealed that the ceRNA pairs displayed significantly increased competition in samples with low miRNA expression (P=1.108×10⁻¹³; Fig. 2B). This is understandable, as mRNAs will be fully suppressed when the interacting miRNA is abundantly expressed; however, when the miRNA molecules are relatively sparse, the targeted mRNAs will only be partially suppressed, and the upregulation of one mRNA will further weaken the miRNA suppression for its partner, resulting in strong co-expression of the ceRNA pairs. Subsequently, the present study analyzed the modular structure of the OSceNet. A total of 13 modules were identified using cFinder software (22). Analysis of the genes in the modules showed that 72% of the nodes from the OSceNet were part of at least one module, thus indicating that the ceRNA pairs tend to interact with other pairs, as a module, in osteosarcoma (Fig. 2C).

To identify the modules with the most robust ability to predict recurrence-free survival for patients with osteosarcoma, the present study randomly divided the 37 samples into two groups: The training dataset, which contained 19 samples, and the test dataset, which consisted of 18 samples. To avoid the prognostic noise of clinical factors, the division was balanced to obtain two datasets with similar clinical information (Table I). For each of the 13 modules identified, the associated recurrence-free survival prognosis was evaluated using a module-specific risk score model (23,24). The risk score model was constructed according to the linear combination of the expression values of genes in the module weighted by the regression coefficients derived from univariate Cox regression analysis (25,26). For example, module 1 included two risk-associated genes, which had positive regression coefficients, and seven protective genes, which had negative regression coefficients, as identified by univariate Cox regression analysis, and the risk score model for module 1 was as follows: Risk score=[0.237× expression value of synaptonemal complex protein 2-like (SYCP2L)]+[0.215× expression value of cytochrome P450 family 3 subfamily A member 5 (CYP3A5)]+[−0.11× expression value of smoothelin-like 2 (SMTNL2)]+[−0.143× expression value of olfactory receptor family 4 subfamily D member 2 (OR4D2)]+[−0.287× expression value of G antigen 12E (GAGE12E)]+[−0.341× expression value of BTB domain-containing 9 (BTBD9)]+[−0.627× expression value of spectrin β, erythrocytic (SPTB)]+[−2.445× expression value of tumor protein p53-regulated apoptosis-inducing protein 1 (TP53AIP1)]+[−2.485× expression value of proteoglycan 3, pro eosinophil major basic protein 2 (PRG3)]. The present study then calculated the risk score of each module for each sample in the training dataset, and divided the samples into a low-risk group (e.g. n=9 in module 1) and a high-risk group (e.g. n=10 in module 1) with the median risk score (e.g. risk score=−49.048 in module 1) as the cutoff. The prognosis of each module was then validated in the test dataset and complete datasets, using the same procedure with an unchanged risk score model and cutoff. The two modules that most robustly predicted the recurrence-free survival prognosis were then selected, as shown in Figs. 3 and 4.

Prognostic module 1 consisted of two risk genes (SYCP2L and CYP3A5) and seven protective genes (SMTNL2, OR4D2, GAGE12E, BTBD9, SPTB, TP53AIP1 and PRG3) (Fig. 3A and B). Of these, previous studies also revealed CYP3A5 to be a risk gene in osteosarcoma. For example, Dhaini et al reported that the expression levels of CYP3A5 were significantly increased in primary biopsies from patients who developed distant metastatic disease compared with biopsies from patients with nonmetastatic disease (27). Furthermore, Trujillo-Paolillo et al also reported that CYP genes serve an important role in osteosarcoma tumorigenesis at primary and metastatic sites, as well as in treatment response (28). The present study revealed a novel potential mechanism by which CYP3A5 may promote osteosarcoma tumorigenesis via its ceRNA activity.

Compared with the low-risk samples in the training, test and complete datasets, the high-risk samples, as defined by module 1 exhibited a higher rate of recurrence (training, 33.3 vs. 90%; test, 0 vs. 50%; total, 20 vs. 68.2%) and a significantly reduced recurrence-free survival time (log-rank P=0.013, 0.05 and 0.008, respectively; Fig. 3C–E).

Prognostic module 2, as shown in Fig. 4A, contained 11 genes, including six risk genes [solute carrier family 7 member 13 (SLC7A13), XK-related 4 (XKR4), synaptotagmin 1, HFM1, ATP-dependent DNA helicase homolog (HFM1), membrane spanning 4-domains A8 and RNA-binding motif protein, Y-linked, family 1, member J] and five protective genes [solute carrier family 22 member 2, usherin (USH2A), DPY30 domain-containing 1 (DYDC1), sorbin and SH3 domain-containing 2 and low-density lipoprotein (LDL) receptor-related protein 1B (LRP1B)]. From this module, LRP1B, which is a member of the LDL receptor family, is particularly notable, as another member of the LDL receptor family (LRP1/LRP1A) is significant in osteoblast behavior (29). In addition, the present study revealed that LRP1B is a protective gene in osteosarcoma (Fig. 4B). The high- and low-risk samples defined by module 2 also exhibited significantly different associated recurrence-free survival prognosis in the training, test and complete datasets (log-rank P=0.047, 0.026 and 0.012, respectively) (Fig. 4C–E).

To determine whether the association of the two modules with recurrence-free prognosis was independent and robust, a multivariate Cox regression analysis was conducted, including risk score defined by the modules, age, sex, primary disease site and chemosensitivity as covariates. The results indicated that risk score was independently significant and the most valuable prognostic factor (HR=9.03, P=0.0045 for module 1; HR=13.42, P=0.0014 for module 2; Table II). In addition, chemosensitivity of the samples was also associated with recurrence-free prognosis (P=0.03 for module 1, P=0.0068 for module 2; Table II). Therefore, the present study subsequently inspected the ability of the risk score to predict recurrence-free prognosis in the samples stratified by chemosensitivity into favorable (n=14) and unfavorable (n=23) strata. The patients within each chemosensitivity stratum were then classified into high and low risk scores based on modules 1 and 2, respectively. For the unfavorable samples, the results demonstrated that the risk score variables derived from modules 1 or 2 had significant prognostic value for recurrence-free survival (log-rank P=1.99×10⁻³ and 5.05×10⁻⁴, respectively, Fig. 5). For the favorable chemosensitivity stratum, the prognostic impact of the modules was reduced (Fig. 5). These results indicated that the prognostic power of the module-based risk score was also chemosensitivity-independent, particularly for patients with unfavorable chemosensitivity. For patients with unfavorable chemosensitivity, early diagnosis may be possible using the risk score. In addition, combining chemosensitivity and the risk score may provide more accurate prognostic information for patients.

The present study revealed that the two modules were independent prognostic factors for recurrence-free survival in osteosarcoma. No common genes existed common between the two modules. When considering upstream miRNAs that regulate the two modules, the same miRNA, hsa-miR-335-3p, was shared by all nodes in modules 1 and 2 (Fig. 6A). Using the median relative expression value for hsa-miR-335-3p across all samples, the samples were divided into high and low hsa-miR-335-3p expression groups. In the two modules, the co-expression coefficients for ceRNA pairs in the low group were significantly higher than those in the high group (P=0.002, <2.2×10⁻¹⁶, respectively), thus suggesting that the two modules are highly regulated by hsa-miR-335-3p (Fig. 6B).

When analyzing the genomic location of the ceRNAs in the modules, it was revealed that the majority of the ceRNA pairs regulated each other across chromosomes, as only three out of 21 and two out of 45 ceRNA pairs were located on the same chromosome for module 1 and module 2 (SYCP2L and BTBD9 on chromosome 6, TP53AIP1 and PRG3 on chromosome 11, OR4D2 and SMTNL2 on chromosome 17 for module 1; HFM1 and USH2A on chromosome 1, XKR4 and SLC7A13 on chromosome 8 for module 2), respectively, which is consistent with previous studies (21,30). In addition, the miRNA-mRNA regulation had a similar result, as only one miRNA-mRNA pair (hsa-miR-335-5p and CYP3A5) was located on the same chromosome in module 1 (Fig. 6C). By combining the expression levels of the genes from the two modules and hsa-miR-335-5p, a novel prognostic model for recurrence-free survival was generated with enhanced power. In this novel risk score model, samples in the high-risk group had a significantly reduced recurrence-free survival time compared with those in the low-risk group (median recurrence-free survival, 28.9 vs. 151.0 months, log-rank P=0.008; Fig. 6D).

The prognostic value for modules 1 and 2, and the combination of the modules with hsa-miR-335-5p, was further validated using ROC curves. The patients with and without recurrence were used as true positive and true negative data-sets, respectively. Good area under the curve (AUC) values were obtained (Fig. 7A). To avoid model overfitting to the training and test datasets, the group division process was repeated with 1,000 iterations, and the AUC values for the three models were evaluated in each training and test dataset. The median AUC values ranged between 0.60 and 0.71 across the different dataset divisions (Fig. 7B). These results indicated that the modules may possess a robust ability to predict the recurrence-free survival of patients with osteosarcoma.