Gene expression datasets were retrieved from Gene Expression Omnibus (GEO) using the key words 'osteosarcoma', 'metastasis', and 'Homo sapiens' by the end of April 28th, 2016. Five datasets (Table I) that met the following criteria were collected: i) gene expression data; ii) OS; and iii) information about metastasis was described.

The gene expression datasets GSE14359 and GSE14827 were acquired with Affymetrix platforms. Thus, background correction and normalization were performed with package 'affy' of R. Missing values were filled with the median method. Background correction was carried out using the MSA method. Normalization was achieved with the quantiles method.

For the other 3 gene expression datasets, probes were then mapped into genes. Probes corresponding to the same gene were averaged as the final expression value of the gene. Normalization was performed with package limma (18) of R.

Prior to meta-analysis, the 5 gene expression datasets were assessed using package MetaQC (19), as well as principal component analysis (PCA) and standardized mean rank (SMR). The following criteria were applied in the assessment: i) internal quality control (IQC), homogeneity test of gene expression profiles among datasets; ii) external quality control (EQC), homogeneity test of gene expression profiles with pathway database; iii) accuracy quality control (AQC), accuracy of featured genes (AQCg) or pathways (AQCp); and iv) consistency quality control (CQC), consistency in ranking of featured genes (CQCg) or pathways (CQCp).

DEGs were screened out using the MetaDE.ES from package MetaDE. This method first tested the heterogeneity of gene expression value in various platforms with three statistic parameters: τ², Q-value and Qpval. It then tested differential expression of genes between different groups with statistic parameters P-value. To ensure the homogeneity of featured genes, τ²=0, Qpval>0.05 and P<0.05 were set as the cut-offs.

The PPI information was downloaded from the database of protein, chemical, and genetic interactions (BioGRID, http://thebiogrid.org/) (20), Human Protein Reference Database (HPRD, http://www.hprd.org/) (21) and Database of Interacting Proteins (DIP, http://dip.doe-mbi.ucla.edu/dip/Main.cgi) (22). The protein products of the DEGs were mapped into the whole network and the PPI network for the DEGs was then acquired. The network was then visualized with Cytoscape (23).

The featured genes were screened out from the DEGs using BC, which reflected hubness of the node in the PPI network. The BC was calculated as follows: where σ_st is the shortest path from s to t; σ_st(ν) is the number of the shortest path from s to t passing through node V; BC score is between 0 and 1, and the greater BC score indicates the higher degree of hubness.

The dataset GSE21257 was selected as the training set. Genes were ranked based upon the BC value and top 10 genes were selected out to train the SVM classifier. An increment of 10 genes were added into the classifier until the metastatic OS could be totally separated from non-metastatic OS. These DEGs were regarded as featured genes. Two-way clustering of sample and gene expression was applied on the featured genes and the result was visualized with a heatmap. The other 3 gene expression datasets were used as the validation set. Sensitivity (Se), specificity (Sp), positive predictive value (PPV), negative predictive value (NPV) and the area under the curve (AUC) were calculated to evaluate the SVM classifier.

Gene Ontology (GO) biological pathways and Kyoto Encyclopedia of Genes and Genomes (KEGG) pathways related to featured genes were identified using runHyperGO and runHyperKEGG from package EMA of R. The significance was calculated using Fisher's exact test as follows: where N is the total number of gene; M is the number of gene in the pathway; K is the number of featured gene.

Quality control test was applied on the 5 gene expression datasets and the results are listed in Table II. As shown in the PCA result (Fig. 1), the dataset GSE9508 far deviated from the other 4 datasets, which was assessed as poor quality. Thus, the dataset GSE9508 was excluded from the following analysis.

According to the criteria described above, a total of 353 DEGs were identified from the 4 datasets. The top 10 DEGs (ranked by P-value) are listed in Table III, including chromosome 13 open reading frame 18 (C13orf18), activating transcription factor 3 (ATF3), cadherin 5 (CDH5), epoxide hydrolase 2 (EPHX2), cadherin 3 (CDH3), ectonucleoside triphosphate diphosphohydrolase 5 (ENTPD5), forkhead box A1 (FOXA1), ephrin B2 (EFNB2), bromo adjacent homology domain containing 1 (BAHD1) and chromosome 10 open reading frame 76 (C10orf76).

Non-DEGs interacting with no less than 5 DEGs were also included in the PPI network. Finally, a PPI network containing 164 nodes and 272 edges was obtained (Fig. 2). The distribution of degree is shown in Fig. 3. Most genes (168 genes) showed a small degree (Log transformed degree <1), while only 2 genes had a Log transformed degree of >4. Therefore, similar to most biological networks, this PPI network exhibited scale-free property. These genes with a high degree were hub genes and may thus play important roles in the development of diseases.

BC was calculated for each node and the top 10 genes are listed in Table IV. These were estrogen receptor 1 (ESR1), E1A binding protein p300 (EP300), tumor protein p53 (TP53), SMAD family member 3 (SMAD3), epidermal growth factor receptor (EGFR), fibronectin 1 (FN1), SRC proto-oncogene, non-receptor tyrosine kinase (SRC), protein kinase Cα (PRKCA), DEAD-box helicase 5 (DDX5) and erythropoietin receptor (EPOR).

Top genes were selected out to train the SVM classifier with the dataset GSE21257. As shown in Fig. 4A, the accuracy increased from 75 to 100%, while more genes were included in the classifier. The accuracy reached 100%, while the top 100 genes (containing 64 DEGs; data not shown) were included in the classifier. Therefore, the top 64 featured genes (data not shown) were selected to construct the SVM classifier. The classifier could separate metastatic OS from non-metastatic OS in the dataset GSE21257 (Fig. 4B).

Hierarchical clustering was performed for the 53 samples from the training dataset GSE21257 using the 64 featured genes and the results are shown in Fig. 5. The 64 featured genes could well distinguish metastatic OS samples from non-metastatic OS samples.

The SVM classifier was validated using the other 3 independent datasets GSE32981, GSE14827 and GSE14359. The accuracies were high in all datasets, with a prediction accuracy of 100% in dataset GSE32981 (Fig. 6A) and GSE14359 (Fig. 6C) and 92.6% in dataset GSE14827 (Fig. 6B). Two non-metastatic OS samples were predicted wrong as metastatic samples in dataset GSE14827. The accuracy, sensitivity, specificity, PPV, NPV and AUC are listed in Table V.

A total of 11 GO biological pathways (Fig. 7) and 5 KEGG pathways (Table VI) were identified in the 64 featured genes. Some pathways were significantly over-represented in the featured genes, such as regulation of cell proliferation, regulation of apoptosis, pathways in cancer, regulation of actin cytoskeleton and the TGF-β signaling pathway.