Bioinformatics analysis of differentially expressed pathways related to the metastatic characteristics of osteosarcoma.

In this study, gene expression data of osteosarcoma (OSA) were analyzed to identify metastasis-related biological pathways. Four gene expression data sets (GSE21257, GSE9508, GSE49003 and GSE66673) were downloaded from Gene Expression Omnibus (GEO). An analysis of differentially expressed genes (DEGs) was performed using the Significance Analysis of Microarray (SAM) method. Gene expression levels were converted into scores of pathways by the Functional Analysis of Individual Microarray Expression (FAIME) algorithm and the differentially expressed pathways (DEPs) were then disclosed by a t-test. The distinguishing and prediction ability of the DEPs for metastatic and non-metastatic OSA was further confirmed using the principal component analysis (PCA) method and 3 gene expression data sets (GSE9508, GSE49003 and GSE66673) based on the support vector machines (SVM) model. A total of 616 downregulated and 681 upregulated genes were identified in the data set, GSE21257. The DEGs could not be used to distinguish metastatic OSA from non-metastatic OSA, as shown by PCA. Thus, an analysis of DEPs was further performed, resulting in 14 DEPs, such as NRAS signaling, Toll-like receptor (TLR) signaling, matrix metalloproteinase (MMP) regulation of cytokines and tumor necrosis factor receptor-associated factor (TRAF)-mediated interferon regulatory factor 7 (IRF7) activation. Cluster analysis indicated that these pathways could be used to distinguish between metastatic OSA from non-metastatic OSA. The prediction accuracy was 91, 66.7 and 87.5% for the data sets, GSE9508, GSE49003 and GSE66673, respectively. The results of PCA further validated that the DEPs could be used to distinguish metastatic OSA from non-metastatic OSA. On the whole, several DEPs were identified in metastatic OSA compared with non-metastatic OSA. Further studies on these pathways and relevant genes may help to enhance our understanding of the molecular mechanisms underlying metastasis and may thus aid in the development of novel therapies.

Introduction

Osteosarcoma (OSA) is the most common primary bone tumor in children and young adults, with 7,104 cases identified during 1999–2008 in the US (1). Although there are a variety of chemotherapeutic and radiation treatments, the survival rate has not greatly improved over the past two decades. Approximately 50% of patients undergo fatal lung metastasis at the late stage of the disease (2,3). The early detection of metastasis is rather critical considering there is a great difference in the survival rate between patients with metastatic OSA (10–20%) and non-metastatic OSA (50–78%) (4,5).

Some achievements have been made in unveiling the molecular mechanisms underlying metastasis. For example, ezrin, an actin-cytoskeleton linker protein, is regarded as a metastatic determinant (6); a high expression of ezrin is associated with a poor outcome in patients with OSA (7). Further studies indicated that ezrin may be involved in metastasis by linking F-actin to the cell membrane following phosphorylation at T567 by protein kinase C) (PKC) (8,9). In addition, heparanase (10), phosphatase and tensin homolog (PTEN) (11) and neuro-pilin2 (12) have also been implicated in the metastasis of OSA. Furthermore, microarray technology has also been widely used for the screening a serial of metastasis-related genes in OSA simultaneously by differential analysis (13,14). However, gene expression signatures show a high variability from different data sets and different literatures, which make them difficult to be used to distinguish metastatic from non-metastatic OSA in the clinic.

Considering that proteins cooperate to exert certain biological functions (15), functional term-based analysis may better distinguish metastatic from non-metastatic OSA. In the present study, gene expression levels were converted into scores of biological pathways. Differentially expressed pathways (DEPs) were identified between metastatic and non-metastatic OSA. The ability of the DEPs to distinguish metastatic OSA from non-metastatic OSA was examined using different gene expression data sets.

Data collection methods

Gene expression data

Gene expression data of OSA were downloaded from Gene Expression Omnibus (GEO, http://www.ncbi.nlm.nih.gov/geo) (16). The following criteria were applied to screen out appropriate gene expression data: i) both metastatic and non-metastatic OSA samples were included; ii) >10 samples were included; iii) human samples.

Four data sets were finally included. Data set GSE21257 (17) contained 34 metastatic OSA tissue samples (from 34 patients, 18.1±12.5 years, 9 females and 25 males) and 19 non-metastatic OSA tissue samples (from 19 patients, 19.4±12.2 years, 10 females and 9 males) that were collected with Illumina Human-6 v2.0 Expression BeadChip (Illumina Inc., San Diego, CA, USA). Data set GSE9508 (18) included 21 metastatic OSA tissue samples (from 11 metastatic patients, 10 of which are analyzed in duplicate, 26.2±19.6 years, 8 females and 3 males), 13 non-metastatic OSA tissue samples (from 7 non-metastatic patients, 6 of which are analyzed in duplicate, 16.6±1.9 years, 2 females and 4 males) and 5 normal controls (66.0±12.9 years, 2 females and 3 males) that were collected with Agilent-Whole Human Genome Oligo Microarray G4112A condensed (Agilent Technologies Inc., Santa Clara, CA, USA). Data set GSE49003 contained 6 metastatic OSA cell samples (KHOS and KRIB, each 3) and 6 non-metastatic OSA cell samples (HOS and U2OS, each 3) that were collected by Illumina Human HT-12 V3.0 Expression BeadChip (Illumina Inc.). Data set GSE66673 contained 12 metastatic OSA cell samples [LM5 cells and Ki-ras transformed HOS cells (143B), each 6] and 12 non-metastatic OSA cell samples(SAOS-2 and HOS, each 6) that were collected by Agilent-028004 SurePrint G3 Human GE 8×60K (Agilent Technologies Inc.).

Data pre-processing and differential analysis of GSE21257

According to annotation files, probes were mapped into genes. Probes mapping to the same gene were averaged as the final expression level of the gene (19). Genes with >20% missing values were excluded and others were filled with the average value, as previously described (20).

Following normalization and log2 transformation, differential analysis was performed using the significance analysis of microarray (SAM) method based on the t-test, as previously described (21). This method can reduce the false-positive by adjusting the p-value to the false discovery rate (FDR) by multiple testing. The relative difference d was calculated as follows: where d is the relative difference, x1′ and x2′ are average gene expression levels at two different statuses, and s is the variance. Genes with a fold change in expression of >2 and a FDR of <0.05 were considered as the differentially expressed genes (DEGs).

Functional profile establishment

Pathway information was downloaded from MsigDB (www.broad.mit.edu/gsea/msigdb/) (22). Gene expression levels were converted into scores of functional terms (such as pathways) using the Functional Analysis of Individual Microarray Expression (FAIME) algorithm (23). First, gene expression levels were converted into rank-related weights. High expression levels meant high weights: where g and s represent gene g and sample s, respectively; r represents the rank of gene g in sample s in ascending order, and |G| represents the total number of genes.

Normalized centroid (NC) is defined as the average weight of a gene set. For a certain functional term Gene Ontology (GO), it has a gene set GOi. For a certain gene g,g∈GOi. G/GOi represents the complementary set of GOi. The NCs for GOi and G/GOi were calculated as follows:

The score of the functional term GO in sample s was defined as the difference between the two NCs as follows: where |GOi| represents the total number of genes in the functional term i. G/GOi is the complementary set of GOi. is the score of functional term i in sample s. FPs represents the scores of all functional terms in sample s. All samples (FP, FP … FP) constitute the functional expression matrix FP.

After the calculation of the functional expression matrix, the t-test was used to screen the DEPs between metastatic and non-metastatic OSA. The p-value was corrected by the Bonferroni method and the adjusted p<0.05 was set as the threshold value.

Prediction of metastasis

A support vector machines (SVM) model was trained with differentially expressed functional terms using function the svm of package e1071 in R (24). The functional expression matrices of data sets GSE9508, GSE49003 and GSE66673 were calculated and the model was then used to predict the metastasis of OSA samples in the 3 data sets. The accuracy (ACC) was calculated as follows: where 'TP' (true-positive) represents right prediction, while 'FP' (false-positive) represents wrong prediction.

Principal component analysis (PCA) was performed to confirm the differential functions of the DEGs and DEPs for the metastatic and non-metastatic OSA samples.

Results

A total of 24,999 genes was identified in data set GSE21275, 3,674 genes in data set GSE9508, 25,159 genes in data set GSE49003 and 18,539 genes in data set GSE66673. The box plots of these 4 data sets are presented in Fig. 1. A good performance of normalization was achieved.

Figure 1

Box plots of gene expression data from data sets (A) GSE21257 and (B) GSE9508. A good performance of normalization was achieved. Box plots of gene expression data from data sets (C) GSE49003 and (D) GSE66673 (D). A good performance of normalization was achieved.

Identification and biological functions of DEGs

Differential analysis was performed between the metastatic and non-metastatic OSA samples using the SAM method. According to the criteria (fold change >2 and FDR <0.05), a total of 616 down-regulated and 681 upregulated genes were identified in data set GSE21257. However, the results of PCA indicated that the DEGs could not be used to not distinguish the metastatic from non-metastatic OSA samples (Fig. 2), suggesting that other classification features were required.

Figure 2

Results of principal component analysis (PCA) based upon the 1,297 differentially expressed genes (DEGs). Red circles represent metastatic osteosarcoma (OSA) samples, while blue circles represent non-metastatic OSA samples.

Functional expression matrix and identification of DEGs

A total of 4,722 pathways (i.e., gene sets) were obtained from MsigDB and 4,721 gene sets which had expression level information for >3 genes were retained. Scores of pathways were calculated for each sample from data set GSE21257 and then a functional expression matrix was obtained. As shown in Fig. 3, a good normalization of pathway scores was achieved, which eliminated the difference in gene expression levels.

Figure 3

Distribution of F score in the 53 samples of the data set GSE21257.

After t-test analysis and Bonferroni correction, a total of 14 DEPs was revealed (adjusted p-value <0.05) (Table I), such as AML1-MTG8 fusion gene, NRAS signaling, Toll-like receptor (TLR) signaling, matrix metalloproteinase (MMP) regulation of cytokines and tumor necrosis factor receptor-associated factor (TRAF)-mediated interferon regulatory factor 7 (IRF7) activation.

Table I

Differentially expressed pathways between metastatic and non-metastatic osteosarcoma.

Pathway term	p-value	Adjusted p-value
Targets of AML1-MTG8 fusion	1.18E-06	0.006
Leishmania infection	3.09E-06	0.015
NRAS signaling	3.46E-06	0.016
T lymphocyte and NK progenitor	3.85E-06	0.018
Endogenous pathway	4.63E-06	0.022
Liver cancer metastasis	4.69E-06	0.022
Trafficking and processing of endosomal TLR	5.94E-06	0.028
Chronic lymphocytic leukemia	6.18E-06	0.029
TARF6-mediated IRF7 activation in TLR7/8 or 9 signaling	7.23E-06	0.034
Thyroid cancer cluster 4	7.31E-06	0.035
Tretinoin response	8.19E-06	0.039
NRAS vs. stromal stimulation	8.61E-06	0.041
Aging cerebellum	8.64E-06	0.041
MMP regulation of cytokines	1.01E-05	0.048

NK, natural killer; TLR, Toll-like receptor; IRF7, interferon regulatory factor 7; MMP, matrix metalloproteinase.

These DEPs contained 400 genes. Compared with the 1,297 DEGs, there were 93 overlapping genes. None of the 14 pathways was significantly over-represented in the 93 DEGs. This finding thus suggested that pathway-based differential analysis was more robust than gene-based analysis at unveiling the molecular mechanisms underlying the metastasis of OSA. As expected, and as shown in Fig. 4, bidirectional hierarchical clustering implied that the metastatic OSA samples were well separated from the non-metastatic OSA samples by these 14 DEPs.

Figure 4

Results of cluster analysis with the 14 differentially expressed functional terms.

Prediction accuracy

An SVM model was trained with the 14 DEPs. The functional expression matrices of data set GSE9508, GSE49003 and GSE66673 were calculated, and the model was then used to predict the metastasis of OSA samples from these 3 data sets. The accuracy for data set GSE9508 was 91%, that for data set GSE49003 was 66.7% and the accuracy for GSE66673 was 87.5%. The detailed information is listed in Table II.

Table II

Detailed information of the prediction results.

Data set	Metastatic samples	Non-metastatic samples	TP/FP for metastatic samples	TP/FP for non-metastatic samples	Accuracy (%)
GSE9508	21	13	21/0	10/3	91
GSE49003	6	6	4/2	4/2	66.7
GSE66673	12	12	10/2	11/1	87.5

TP, true-positive; FP, false-positive.

The results of PCA for the 3 data sets are shown in Fig. 5. A good prediction was also achieved with the 14 DEPs, further suggesting the reliability of the DEPs and our method.

Figure 5

Results of principal component analysis (PCA) based upon the 14 differentially expressed functional terms. (A) GSE9508, (B) GSE49003 and (C) GSE66673. Metastatic osteosarcoma samples are shown in red while non-metastatic osteosarcoma samples are shown in blue.

Discussion

Gene-based differential analysis of gene expression data is a useful method to identify disease-related genes. However, numerous genes are identified in a study and it takes time to further screen out critical genes. Proteins interact with each other to fulfill certain biological functions. Therefore, pathway-based differential analysis could disclose pathways disturbed in disease. In present study, gene expression levels were converted into pathway scores and 14 DEPs were identified in metastatic OSA, such as NRAS signaling, TLR signaling, MMP regulation of cytokines and TRAF-mediated IRF7 activation. Their ability to distinguish metastatic OSA from non-metastatic OSA was proven with 3 different gene expression data sets.

The expression of most TLRs is heterogeneous in OSA cell lines, apart from TLR-2 and TLR-4 (25). The expression levels of TLRs are differentially influenced by p53 and DNA stresses in cancers (26). That is why a pathway-based differential analysis is more robust. A previous study indicated that the expression levels of TLR3, TLR4 and TLR9 have clinical interest as indicators of tumor aggressiveness in breast cancer (27). TLRs also mediate metastasis induced by pancreatic adenocarcinoma upregulated factor (28). The functional expression of TLR9 is related to the metastatic potential of human lung cancer cells (29). Therefore, we speculated that TLR signaling plays a similar role in OSA metastasis and some members of TLRs (such as TLR-7, -8 and -9) may play critical roles in the metastasis.

Several members of TRAFs have been shown to be involved in the growth and invasion of OSA. TRAF6 regulates the proliferation, apoptosis and invasion of OSA cells (30). The knockdown of TRAF4 expression suppresses OSA cell growth in vitro and in vivo (31). The study by Yao et al indicated that TRAF4 can also enhance OSA cell invasion through the AKT signaling pathway (32). IRF7 pathways have also been shown to be involved in metastasis. Bidwell et al reported that the silencing of IRF7 in breast cancer cells promoted bone metastasis through immune escape (33).

MMPs also participate in OSA invasion (34). The effect of S100A4 on metastasis is mediated by MMPs (35). Previous studies have confirmed that MMP-2 and MMP-9 play critical roles in the modulation of metastasis in OSA (36,37). The neuroblastoma RAS viral oncogene homolog (NRAS) mutation status is an independent prognostic factor in metastatic melanoma (38). The loss of NRAS induces the metastatic conversion of Rb1-deficient neuroendocrine thyroid tumors (39). We hypothesized that NRAS signaling plays a similar role in the metastasis of OSA.

Furthermore, immune-related pathways were also differentially expressed in metastatic OSA, including the T lymphocyte and natural killer (NK) progenitor and IRF7 pathway. Immune cells regulate metastatic progression in solid tumors (40). T cell based immune surveillance plays an important role in preventing OSA metastasis (41). NK cells are also involved in metastasis (42,43). The modulation of NK cells may be a method for treating OSA (44).

In conclusion, several DEPs were identified in metastatic OSA. Limited to the disclosed gene expression data sets of OSA, these results warrant further experimental confirmation. However, the accuracy for the prediction of metastasis suggested that the identified pathways are valuable clues. Further studies on these pathways would enhance our understanding of the mechanisms underlying the metastasis of OSA. Potential biomarkers and therapeutic targets may also be identified from these pathways.