The GSE14359 (8) gene expression profile data were acquired from the Gene Expression Omnibus (http://www.ncbi.nlm.nih.gov/geo/), which was based on the platform of the GPL96 [HG-U133A] Affymetrix Human Genome U133A Array. This dataset contains 10 conventional OS samples from the femur or tibia (two replicates each) from five consenting patients with grade 2–3 conventional OS between 7 and 74 years of age; eight OS lung metastasis tumour samples (two replicates each) from four consenting patients with a grade 1–3 OS lung metastatic tumour; and two non-neoplastic primary osteoblast cell samples with limited life span in vitro from one patient (two replicates). These 10 conventional OS samples and two non-neoplastic primary osteoblast samples were selected for further analysis.

The CEL files and probe annotation files were downloaded, and the gene expression data of all samples were preprocessed via background correction, quantile normalization and probe summarization using the Gene Chip Robust Multi Array algorithm (12) in the Affy software package (version 1.32.0; http://www.bioconductor.org/packages/release/bioc/html/affy.html) (13).

The Linear Models for Microarray Data package (version 3.10.3; http://www.bioconductor.org/packages/2.9/bioc/html/limma.html) (14) of R was used to identify genes, which were significantly differentially expressed in the conventional OS samples. The raw P-value was adjusted by the Benjamin and Hochberg method (15), and only genes meeting the cut-off criteria of |log₂ fold-change|>1 and adjusted P<0.01 were selected as DEGs.

Hierarchical clustering is a common method used to determine clusters of similar data points in multidimensional spaces (16). The pheatmap package (version 1.08; https://cran.r-project.org/web/packages/pheatmap/) (17) was used to perform hierarchical clustering via joint between-within distances for the DEGs in the conventional OS and non-neoplastic primary osteoblasts samples.

The Database for Annotation, Visualization and Integrated Discovery (DAVID) provides a set of comprehensive functional annotation tools, which can be used to identify the biological meanings of abundant genes (18). P<0.01 was used as the cut-off criterion for GO and Kyoto Encyclopedia of Genes and Genomes pathway enrichment analysis using DAVID (version 6.7; https://david-d.ncifcrf.gov/), based on the hypergeometric distribution algorithm.

The Search Tool for the Retrieval of Interacting Genes database (version 10.0; http://string-db.org/), which provides experimental and predicted interaction information (19), was used to analyse the PPIs for DEGs by calculating the combined score, and a score >0.4 was selected as the cut-off criterion. Subsequently, the PPI network of the upregulated and downregulated DEGs was visualized using Cytoscape (version 3.2.0; http://cytoscape.org/) (20).

The network modules were obtained based on Molecular Complex Detection (MCODE) analysis (21) of the original PPI networks. The default parameters (degree cut-off, 2; node score cut-off, 0.2; K-core, 2) were used as the cut-off criteria for the network module screening.

In order to obtain further information on the gene functions and identify pathways closely associated with the DEGs, functional annotation analysis and subsequent pathway enrichment analysis of the network module with the highest MCODE score were performed using DAVID, with a P<0.01 cut-off.

The potential miRNAs for upregulated DEGs in the most significant pathway were predicted using the miRDB database (version 1.24.0; http://www.bioconductor.org/packages/2.8/bioc/html/maDB.html) (22), with a cut-off for the target score of ≥60. The binding sites of miRNAs in the human mRNAs > 800 were abandoned. The integrated miRNA-DEG regulatory network was then visualised with Cytoscape.

Following the data preprocessing, 11,107 probes were obtained. Based on the cut-off criteria, a total of 987 DEGs were screened from the conventional OS samples, including 317 upregulated genes and 670 downregulated genes. The hierarchical cluster analysis of the data revealed that it was possible to use the DEGs to accurately distinguish the conventional OS samples from the non-neoplastic primary osteoblast cell samples (Fig. 1).

According to the GO functional annotation, the upregulated DEGs were predominantly enriched in GO terms associated with DNA replication, including MCM3, replication factor C (RFC)5, replication protein A3 (RPA3) and flap endonuclease 1 (FEN1), and cell cycle, including cyclin-dependent kinase 1 (CDK1), NDC80 kinetochore complex (NDC80), BUB1 mitotic checkpoint serine/threonine-protein kinase (BUB1) and mitotic arrest deficient 2 like 1 (MAD2L1). A number of downregulated DEGs, including caveolin 1 (CAV1), cadherin 13 (CDH13), vascular endothelial growth factor C (VEGFC) and transforming growth factor β receptor 3 (TGFBR3), were relevant to blood vessel development, whereas epidermal growth factor receptor (EGFR), TP53, VEGFB and MAPK1 were associated with the regulation of cell proliferation (Table IA).

According to the results of the pathway enrichment analysis, the upregulated DEGs were predominantly enriched in seven pathways. In accordance with the GO term analysis, the DNA replication pathway, including RFC2, RFC3, RFC4 and RFC5, and cell cycle pathway, including CDK1, minichromosome maintenance complex component 3 (MCM3) and BUB1, were also enriched in the upregulated genes. The downregulated DEGs were predominantly enriched in the focal adhesion, including CAV1, collagen type VI α1 (COL6A1), thrombospondin 1 (THBS1) and EGFR, and p53 signalling pathways, including TP53, Fas cell surface death receptor (FAS) and TP53 apoptosis effector (PERP), as shown in Table IB.

The PPI network for the upregulated and downregulated DEGs consisted of 442 pairs of PPIs. The degrees of DEGs, including CDK1, MAD2L1, NDC80, non-SMC condensin I complex subunit G (NCAPG), BUB1, centromere protein F (CENPF) and kinesin family member 11 (KIF11), were >17 (Table II), indicating that they were important genes in OS.

A total of 10 network modules were obtained from the PPI network using the default criteria, and the module with the highest score contained 16 nodes and 102 edges. In this module, CDK1 interacted with other DEGs, including MAD2L1, BUB1, NCAPG, NDC80 and CENPF (Fig. 2).

The functional enrichment analysis of the module with the highest score showed that the majority of the DEGs in this module were predominantly associated with the cell cycle. Certain DEGs, including CDK1, MAD2L1, BUB1 and NDC80, were implicated in mitosis and the M phase of mitotic cell cycle; other DEGs, including Rac GTPase-activating protein 1 (RACGAP1) and MAD2L1, were correlated with cell cycle process (Table IIIA). CDK1, MAD2L1, BUB1 and aurora kinase A (AURKA) were significantly enriched in the oocyte meiosis pathway (Table IIIB).

The miRNA-DEG regulatory network contained 63 miRNAs, nine of their corresponding DEGs and 16 DEGs, which interacted with these nine DEGs. DNA polymerase ζ subunit 3 (POLE3) was regulated by 18 miRNAs, including miR-4310, miR-4680-3p, miR-583 and miR-4269; RFC3 was regulated by 16 miRNAs, including miR-802 and miR-649; RFC3 and RFC5 were modulated by miR-522-3p and miR-224-3p. In addition, RFC2, RFC3, RFC4 and RFC5 interacted with each other (Fig. 3).