Regulatory network of differentially expressed genes in metastatic osteosarcoma

Introduction

Osteosarcomas (OS) are among the most frequently occurring secondary malignancies in childhood cancer (1). OS most often originates in the metaphyses of long bones in adolescents and young adults (2). During the past 30 years, an optimal treatment strategy for OS has been developed, which consists of multi-agent chemotherapy and aggressive surgical resection of all sites of disease involvement (3). However, ~80% of patients with localized OS develop metastatic disease following surgical resection (4). Patients with primary metastatic OS are a heterogeneous group, and a five-year event-free survival rate of up to 75% is reported for patients presenting with unilateral lung metastases (5). Furthermore <20% of patients with high-grade osteosarcoma are clinically diagnosed with metastatic disease at the initial diagnosis, and long-term survival rates of patients with metastatic OS range between 10–40% (6,7). Therefore, due to the high rate of systemic spread, complete recovery following surgical treatment alone is rare (8).

Recently, microarray analysis (9) has been widely used in screening for differentially expressed genes (DEGs), in order to identify potential genes that may be investigated for the treatment of various human diseases (10,11). Previous studies have focused on the microarray analysis of OS. A previous study used genome-wide cDNA microarrays to investigate the transcriptome profile of the human Saos-2 and U-2 OS cell lines. Genes associated with focal adhesion were shown to be differentially regulated in the two cell lines (12). Luo et al (13), identified a total of 1,836 DEGs in OS. In addition, a previous study detected a total of 35 aberrantly expressed genes in three cell lines of OS (IOR/OS9, IOR/OS10, and IOR/OS15), eight were upregulated and 27 were downregulated, as compared with the expression levels in osteoblasts (14). Another DEG analysis yielded 75 upregulated genes and 97 downregulated genes in osteoblastic, as compared with non-osteoblastic OS samples (15). Based on comparative genomic hybridization, chromosomal imbalances have also been shown to exist in OS (16). Therefore, DNA microarray analysis may be considered an effective approach for the identification of genes associated with OS, and may provide potential treatment strategies for OS.

The present study used a microarray analysis to identify the DEGs between metastatic and non-metastatic OS samples, using the BioConductor package in R language. Functional annotation and pathway enrichment analyses of the DEGs were conducted using the MAS system and the Kyoto Encyclopedia of Genes and Genomes (KEGG) pathway, respectively. In addition, the co-expression of DEGs was analyzed using Pearson’s correlation coefficient. Protein-protein interaction (PPI) networks of the co-expressed genes were constructed using Cytoscape and the submodules of the network were selected for by MCODE. The results of the present study may provide increased understanding of the underlying molecular mechanisms of metastatic OS, and identify novel therapeutic targets for its diagnosis and treatment.

Materials and methods

Samples

The gene expression profile GSE9508 (17) was downloaded from the public functional genomics database Gene Expression Omnibus (GEO, http://www.ncbi.nlm.nih.gov/geo/). A total of 23 specimens, which included 11 human metastatic OS samples, seven non-metastatic OS samples and five normal samples, were obtained. The GSE9508 expression profile is based on the GPL6076 Whole Human Genome Oligo Microarray G4112A platform (Agilent Technologies, Inc., Santa Clara, CA, USA).

Data pretreatment and analysis of DEGs

Using the BioConductor package (http://bioconductor.org/) and Gene Spring software (Silicon Genetics, Redwood City, CA, USA) in R language, the probe-level data were converted into expression values. The expression values of all of the probes in each sample were reduced to a single value by determining the average expression value. Missing data were imputed and quantile normalization for complete data was performed, as previously described (18), using the preprocessCore package in R language (19). When numerous probes were mapped to one gene, the mid-value of the data was defined as the expression level of the gene. However, when numerous genes were mapped by one probe, this probe was considered to lack specificity, and was removed from the analysis.

The normal samples were classed as the controls and the normalized data were analyzed using a t-test (20,21) implemented in LIMMA package (22). The differential genes in the metastatic and non-metastatic OS samples were considered to be those with a P-value <0.05. The differential genes that were uniquely expressed in the metastatic OS samples, as compared with the non-metastatic OS samples, were considered to be DEGs.

Functional enrichment and pathway enrichment analyses

Gene Ontology (GO) (23) and KEGG pathway (24) enrichment analyses of DEGs were conducted using the Molecule Annotation System 3.0 (MAS, http://bioinfo.capitalbio.com/mas3/). GO functional enrichment analysis encompasses three categories: Biological process, cellular component and molecular function, with a cutoff criteria of gene count >2 and P<0.01. The KEGG pathway analysis was performed using the same method and parameter settings (count >2 and P<0.01), to identify the pathways that the DEGs are involved in.

Co-expression analysis of the DEGs

The co-expression of DEGs in the metastatic OS samples was also analyzed. Based on the previously obtained DEG, the co-expression of numerous pairs of DEGs in metastatic OS was identified by Pearson’s correlation coefficient (PCC) (5,25,26). A pair of DEGs with PCC ≥0.98 and P<0.01, were identified as significantly co-expressed genes, under specific conditions.

Construction of a co-expression network and identification of functional modules

The human PPI network was downloaded from the Human Protein Reference Database (http://www.hprd.org/) (27). The specific interactional patterns of DEGs in metastatic OS were extracted according to the human PPI network, and the dysfunctional protein networks associated with metastatic OS were constructed by integrating the specific interactional patterns of DEGs and the gene co-expression information. The genes which belonged to the top 10% (nodal points ≥ 6) were screened for further analysis.

The topological characteristics of the co-expression network were examined using Cytoscape MCODE (28), which is a clustering algorithm used for directed or undirected graphs. The MCODE algorithm includes the following steps: Vertex-weighting, complex prediction and optional post-processing (29). The weighting scheme defined a measure of local density for a vertex’s neighborhood. Complexes with a high vertex weight were then used as seed and the complex neighbor vertices were checked to determine whether they were a part of the complex. Post-processing was performed with a minimum degree of 3. The submodules were then obtained from the regulatory network.

Results

Screening DEGs in metastatic OS

Pretreatment and standardization of the data from all of the available samples was conducted, and the DEGs of non-metastatic and metastatic OS samples were compared with the normal samples. A total of 965 characteristic genes were then identified as DEGs in metastatic, as compared with non-metastatic, OS samples (P<0.05).

GO clustering and KEGG pathway analysis of the DEGs

To explore the function of the DEGs in metastatic OS, the DEGs were mapped to the GO database using the MAS system. The results were analyzed based on three categories: Biological process, cellular component and molecular function. The DEGs clustered according to biological process are shown in Table I, and include regulation of transcription, DNA-dependent (P=3.11E-82), signal transduction (P=2.14E-45) and cell adhesion (P=3.65E-32). These results indicate that the majority of DEGs were associated with transcription, signal transduction and cell adhesion. The DEGs clustered according to cellular component were shown to predominantly be associated with the nucleus (P=9.23E-250), cytoplasm (P=2.59E-177) and integral to the membrane (P=4.41E-126) (Table II), thus suggesting that the cellular components in which the DEGs were expressed, was relatively comprehensive. The clustering of DEGs according to molecular function is shown in Table III, and included binding activities, such as protein binding (P=3.43E-232), zinc ion binding (P=2.23E-94) and nucleotide binding (P=2.03E-88).

Table I Gene ontology (GO) functional enrichment of differentially expressed genes according to biological process, in metastatic osteosarcoma.

Table I

Gene ontology (GO) functional enrichment of differentially expressed genes according to biological process, in metastatic osteosarcoma.

GO accession number	Count	Adjusted P-value	Description
GO:0006355	100	3.11E-82	Regulation of transcription, DNA-dependent
GO:0006350	81	3.09E-56	Transcription
GO:0007165	81	2.14E-45	Signal transduction
GO:0055114	37	7.82E-36	Oxidation reduction
GO:0007155	37	3.65E-32	Cell adhesion
GO:0007275	52	3.34E-28	Development
GO:0015031	25	6.24E-18	Protein transport
GO:0044419	24	6.51E-18	Protein amino acid phosphorylation
GO:0006468	18	2.36E-17	Interspecies interaction between organisms
GO:0000122	15	1.51E-16	Negative regulation of transcription from RNA polymerase II promoter

Table II Gene ontology (GO) functional enrichment of differentially expressed genes according to cellular component, in metastatic osteosarcoma.

Table II

Gene ontology (GO) functional enrichment of differentially expressed genes according to cellular component, in metastatic osteosarcoma.

GO accession number	Count	Adjusted P-value	Description
GO:0005634	265	9.23E-250	Nucleus
GO:0005737	236	2.59E-177	Cytoplasm
GO:0016021	174	4.41E-126	Integral to membrane
GO:0016020	191	1.56E-123	Membrane
GO:0005886	125	1.84E-92	Plasma membrane
GO:0005576	90	1.24E-74	Extracellular region
GO:0005829	57	4.75E-55	Cytosol
GO:0005783	50	1.76E-46	Endoplasmic reticulum
GO:0005739	52	3.91E-46	Mitochondrion
GO:0005794	43	1.16E-38	Golgi apparatus

Table III Gene ontology (GO) functional enrichment of differentially expressed genes according to molecular function, in metastatic ostersarcoma.

Table III

Gene ontology (GO) functional enrichment of differentially expressed genes according to molecular function, in metastatic ostersarcoma.

GO accession number	Count	Adjusted P-value	Description
GO:0005515	283	3.43E-232	Protein binding
GO:0008270	109	2.23E-94	Zinc ion binding
GO:0000166	100	2.03E-88	Nucleotide binding
GO:0046872	112	4.23E-72	Metal ion binding
GO:0005524	72	2.55E-66	ATP binding
GO:0003677	70	2.03E-47	DNA binding
GO:0003700	48	5.97E-43	Transcription factor activity
GO:0016740	57	1.43E-41	Transferase activity
GO:0005509	44	3.33E-38	Calcium ion binding
GO:0003723	37	1.10E-33	RNA binding

To further explore the detailed changes to biological pathways in metastatic OS, a KEGG pathway enrichment analysis was conducted to select the altered pathways. The top 10 pathways were identified and are listed in Table IV. The most significantly enriched pathways were associated with the ribosome (P=1.03E-07), axon guidance (P=1.37E-07), cytokine-cytokine receptor interaction (P=3.93E-07) and focal adhesion (P=8.58E-06), which was also enriched in the DEGs clustered according to biological process.

Table IV Kyoto Encyclopaedia of Genes and Genomes (KEGG) pathways of differentially expressed genes in metastatic osteosarcoma.

Table IV

Kyoto Encyclopaedia of Genes and Genomes (KEGG) pathways of differentially expressed genes in metastatic osteosarcoma.

KEGG	Count	Adjusted P-value	Description
KEGG_PATHWAY	10	3.36E-08	Pathogenic Escherichia coli infection-EHEC
KEGG_PATHWAY	10	3.36E-08	Pathogenic Escherichia coli infection-EPEC
KEGG_PATHWAY	12	1.03E-07	Ribosome
KEGG_PATHWAY	13	1.37E-07	Axon guidance
KEGG_PATHWAY	17	3.93E-07	Cytokine-cytokine receptor interaction
KEGG_PATHWAY	8	8.28E-06	Adipocytokine signaling pathway
KEGG_PATHWAY	13	8.58E-06	Focal adhesion
KEGG_PATHWAY	10	2.34E-05	Tight junction
KEGG_PATHWAY	12	5.07E-05	Regulation of actin cytoskeleton
KEGG_PATHWAY	13	5.28E-05	Neuroactive ligand-receptor interaction

PPI network and submodules of DEGs

A gene co-expression analysis was performed in metastatic OS, based on the PCC values. Among the 965 DEGs, a total of 182 gene pairs were identified as co-expressed genes pairs with PCC>0.98 and P<0.01.

In addition, the interaction patterns of DEGs in metastatic OS were obtained by mapping them to the human PPI network. The dysfunctional network associated with metastatic OS was then constructed, by integrating the interaction data of the DEGs and the associations of co-expressed genes (Fig. 1). This network included the interaction of DEGs in metastatic OS samples. In the network, hub nodes included matrix metalloproteinase 1 (MMP1), smoothened (SMO), ewing sarcoma breakpoint region 1 (EWSR1), unnamed protein product 1, sarcolemmal membrane-associated protein (SLMAP), fasciculation and elongation protein ζ-1 (FEZ1).

Figure 1 Protein-protein interaction (PPI) network of coexpressed genes in metastatic osteosarcoma. Green nodes, the differentially expressed genes (DEGs) of the PPI network; yellow nodes, the DEGs only represented within the co-expression network; blue nodes, the DEGs identified in the PPI and gene-coexpression networks.

Four functional submodules, which were dysfunctional in the process of metastatic OS, were then selected using Cytoscape (Fig. 2). MMP1, SMO, brain selective kinase 2 (BRSK2), aldo-keto reductase family 1 member B10 (AKRIB10), as well as other genes were included within these four submodules.

Figure 2 Functional submodules of the protein-protein interaction (PPI) network. Post-processing was performed with a node degree of 3. Green nodes, the differentially expressed genes (DEGs) of the PPI network; blue nodes, the DEGs identified in the PPI and gene-coexpression networks.

Discussion

OS is one of the most frequent secondary cancers that occurs following childhood malignancy, particularly metastatic OS. Due to the low recovery rates and the high incidence of metastatic OS, it is necessary to explore the molecular mechanisms of OS, in order to identify an effective prevention and treatment strategy.

In the present study, 965 DEGs between metastatic and non-metastatic OS samples were identified by comparing their gene expression profiles with those of normal samples. The DEGs were then subjected to functional annotation, based on three categories: Biological process, cellular component and molecular function. A KEGG pathway analysis was also performed. The ribosome and axon guidance were identified as the most significantly altered pathways. Biogenesis and translational control are essential cellular processes associated with ribosomes, which are governed at numerous levels. Numerous tumor suppressors and proto-oncogenes have previously been shown to either affect the formation of the mature ribosome or regulate the activity of translation factors (30). Therefore, the DEGs identified in the present study may have important effects on metastatic OS, through the ribosome pathway.

Axons are guided along specific pathways by attractive and repulsive cues in the extracellular environment (31). The results of the present study suggest that axon guidance may be an important pathway in metastatic OS. Furthermore, the cytokine-cytokine receptor interaction pathway has previously been shown to be present in astrocytomas, and the autoregulation of interleukin-1 and cytokine-receptor interactions was shown to exist in primary human astrocytoma cells (32). In the present study, focal adhesion was identified in the GO biological process analysis of DEGs and the KEGG pathway. Focal adhesion was previously identified as a prominent determinant in cancer initiation, progression and metastasis (33–36).

The present study constructed a PPI network of the co-expressed genes, which was then screened for submodules. Genes, such as MMP1, SMO, EWSR1 and FEZ1, were selected as the hub nodes, suggesting that these genes are associated with the crucial function of the whole network. MMP1 not only existed in the PPI network, but was also identified as one of the co-expressed genes in metastatic OS. The functions of MMP1 have been reported in previous studies. It has been identified as a candidate marker that may be useful for identification of breast lesions, which can develop into cancer (37). MMP1 has also been shown to possess important functions in OS primary tumors and OS metastasis to the lung, which is the predominant site of OS metastasis (38). There has also been a correlation reported between MMP expression and the oncological outcome of OS patients, thus suggesting the prognostic significance of MMPs in OS (39). Furthermore, the regulation of MMP gene expression has vital roles in tumor invasion (40). Hence, MMP1 may have a vital role in metastatic OS.

SMO is a distant relative of the G protein-coupled receptors. It has been shown to mediate the Hedgehog signaling pathway during embryonic development and may initiate ligand independent pathway activation in tumorigenesis (41). The inactivation of SMO has been suggested as a potential target for the treatment of patients with OS (42). EWSR1 represents one of the most commonly involved genes in sarcoma translocations (43), and in a study regarding a fusion transcript in osteogenic sarcoma, it was shown to be closely associated with the molecular mechanisms of small cell osteogenic and Ewing sarcomas (44). EWSR1 gene rearrangements have also been identified in soft tissue myoepithelial tumors (45). FEZ1 is a tumor suppressor gene that maps to chromosome 8p22, which is a frequently deleted chromosomal region in numerous human malignancies (46). A previous study showed that it was associated with p53 (47), it may be suggested that the abnormal expression of FEZ1 may inhibit the normal function of p53, resulting in the occurrence and metastasis of cancer. Furthermore, BRSK2 and AKRIB10 were included in the selected submodules. BRSK2 has been identified as a member of the AMP-activated protein kinase related kinases (48); however, its function in metastatic OS was previously unclear. The results of the present study showed that the molecular function of BRSK2 was associated with nucleotide binding and protein amino acid phosphorylation. Furthermore, research regarding AKRIB10 is limited; however, the present study demonstrated that it was mainly enriched within the cytoplasm, and was involved in the bisphenol A degradation and bile acid biosynthesis pathways. These results may provide novel information regarding the potential mechanisms of metastatic OS.

In conclusion, 965 characteristic DEGs were identified in metastatic, as compared with non-metastatic, OS samples by microarray analysis. Functional annotation and pathway enrichment analyses of the DEGs were also performed, and a total of 182 co-expressed gene pairs were identified in the metastatic OS samples. A regulatory network of co-expressed genes and DEGs was constructed, and the submodules were shown to contain BRSK2, AKRIB10 and other genes. Certain hub nodes identified in the present study, such as MMP1, SMO, EWSR1 and FEZ1, may have the potential to become targets for the diagnosis and treatment of metastatic OS. In addition, BRSK2 and AKRIB10 may have important functions in metastatic OS. The present study provided novel information that may be beneficial for further research regarding the therapy of metastatic OS. However, these results require further confirmation.

Acknowledgements

The present study was supported by the Liaoning Province Science and Technology Plan Project (grant no. 2012408002) and the National Natural Science Foundation (grant no. 81271538).

Abbreviations:

References