Screening of candidate key genes associated with human osteosarcoma using bioinformatics analysis.

The aim of the present study was to identify the key genes associated with osteosarcoma (OS) using a bioinformatics approach. Microarray data (GSE36004) was downloaded from the Gene Expression Omnibus database, including 19 OS cell lines and 6 normal controls. Differentially expressed genes (DEGs) in the OS cell lines were identified using the Limma package, and differentially methylated regions were screened with methyAnalysis in R. Copy number analysis was performed and genes with copy number gains/losses were further screened using DNAcopy and cghMCR packages. Functional enrichment analyses were performed using the Database for Annotation, Visualization and Integrated Discovery online tool, and protein-protein interactions were identified based on information obtained from the Search Tool for the Retrieval of Interacting Genes database. A total of 47 downregulated genes were screened in hyper-methylated regions, including the fragment crystallizable (Fc) region of immunoglobulin E, high affinity I, receptor for; γ polypeptide (FCER1G), leptin (LEP) and feline Gardner-Rasheed sarcoma viral oncogene homolog (FGR). In addition, a total of 17 upregulated genes, including the TPase family, AAA domain containing 2 (ATAD2) and cyclin-dependent kinase 4 (CDK4), exhibited copy number gains, while 5 downregulated genes, including Rho GTPase activating protein 9 (ARHGAP9) and major histocompatibility complex, class II, DO α (HLA-DOA), exhibited copy number losses. These results indicate that hyper-methylation of FCER1G, LEP, and FGR may serve a crucial function in the development of OS. In addition, copy number alterations of these DEGs, including ATAD2, CDK4, ARHGAP9 and HLA-DOA, may also contribute to OS progression. These DEGs may be candidate targets for the diagnosis and treatment of this disease.

Introduction

Osteosarcoma (OS) is the most common type of primary bone malignancy in children and adolescents (1). It is a high-grade malignant tumor with a poor prognosis, and ~20% of patients with OS present with metastases at diagnosis (2,3). Treatment generally involves surgery and adjuvant chemotherapy, and a positive response to chemotherapy is considered to be a prognostic marker for OS, indicating a favorable overall outcome (4,5). Although neoadjuvant chemotherapy has improved the survival rate by a considerable amount, progress in treatment regimens remains at a plateau.

At the molecular level, OS is characterized by gross changes in gene expression and highly heterogeneous karyotypes with variable genomic aberrations (6–8). DNA methylation, particularly at cytosine-phosphate-guanine sites in the promoter region of genes, is a crucial mechanism for the downregulation of gene expression (9,10). Expression levels of Ras association domain family member 1A have been demonstrated to decrease due to promoter methylation in primary OS cell lines (11,12). Hypermethylation of the hypermethylated in cancer 1 promoter has been revealed to be present in 17% of pediatric patients with OS (13). In addition, genetic variations, particularly single nucleotide polymorphisms, may contribute to cancer risk and progression (14). Copy number variations across the whole genome, including deletions, amplifications and duplications, are associated with OS tumorigenesis (15). In addition, the mutation or loss of gene expression of the tumor protein p53 tumor suppressor gene is commonly associated with OS (16,17). Despite these data, there is limited understanding of the molecular pathogenesis of OS, and a lack of good diagnostic and prognostic tools.

A number of previous studies have combined gene expression data with DNA copy number data to screen for tumor-associated genes in OS genetic variants (18,19). The microarray data GSE36004 has been used in the integrative analysis of mRNA expression, DNA methylation and DNA copy number in OS for elucidating dependencies, and the association between genetic and epigenetic alterations in OS (20). In contrast to a previous study (20), the present study aimed to utilize an integrative bioinformatics approach to map the genetic and epigenetic changes in OS, and to identify key genes associated with OS oncogenesis. Microarray data from 19 OS cell lines and normal controls was used to screen differentially expressed genes (DEGs) and differentially methylated regions (DMRs), and to perform copy number analysis. Then, functional enrichment analysis was performed and protein-protein interactions (PPI) were identified to additionally screen for the key genes. The identification of molecular targets that are specific for OS will be critical to the development of novel targeted therapeutic strategies to improve patient outcomes.

Materials and methods

Microarray data

The microarray data GSE36004, which was contributed by Kresse et al (20), were downloaded from the public repository Gene Expression Omnibus (GEO, http://www.ncbi.nlm.nih.gov/geo/) (21). It contained expression profiling data, methylation profiling data and genome variation single nucleotide polymorphism (SNP) profiling data, which were respectively based on Illumina HumanMethylation27 BeadChip (Illumina, Inc., San Diego, CA, USA), Affymetrix Genome-Wide Human SNP 6.0 Array (Affymetrix, Inc., Santa Clara, CA, USA) and Illumina human-6 v2.0 expression BeadChip (Illumina, Inc.) analyses. A total of 25 samples, including 19 OS cell lines and 6 normal controls (osteoblasts and bones), were applied to develop expression profiling data and methylation profiling data, respectively, while only the 19 OS cell lines were utilized to develop SNP profiling data. All the raw data and annotation files were obtained for subsequent analysis.

DEG screening

According to the expression profiling data of GSE36004, DEGs in OS cell lines compared with normal controls were identified using the Limma package (available at http://www.bioconductor.org/packages/release/bioc/html/limma.html) (22) in Bioconductor package version 1.0.2 (23). Significant P-values were adjusted for multiple testing using the Benjamini-Hochberg method (24). A log fold-change (FC) >1 and adjusted P<0.05 were considered to indicate a statistically significant difference. Unsupervised hierarchical clustering of 25 samples was performed using the R package (v2.13.0; R Project for Statistical Computing, Vienna, Austria) (25) based on microarray data.

Analysis of methylated regions

Quantile normalization of the methylation profiling data was first performed using the Lumi package (available at http://www.bioconductor.org/packages/release/bioc/html/lumi.html) (26) in Bioconductor package version 1.0.2, then DMRs were identified using methyAnalysis (27) with minimum P-values adjusted to <0.01.

Copy number analysis

Copy number analysis of the raw data from the SNP profiling data was performed using the crlmm package (http://www.Bioconductor.org/packages/release/bioc/html/crlmm.html) (28) in Bioconductor package version 1.0.2. Then, genes with copy number gains/losses were additionally screened using the DNAcopy (29) and cghMCR package (http://www.bioconductor.org/packages/2.4/bioc/html/cghMCR.html). Genes with copy number gains/losses (0.2 for gain, −0.2 for loss) in >30% samples was the threshold value for the DNAcopy and cghMCR packages.

Functional enrichment analysis of DEGs

Gene Ontology (GO; http://www.geneontology.org) (30) is a tool for the unification of biology functions based on gene annotation information, which primarily consists of biological process (BP), molecular function (MF), and cellular component (CC) analysis. The Kyoto Encyclopedia of Genes and Genomes (KEGG; available at http://www.genome.ad.jp/kegg/) (31) is a pathway-associated database which connects known information on molecular interaction networks. To understand the biological significance of the identified DEGs, the upregulated and downregulated genes were input into the Database for Annotation Visualization and Integrated Discovery (DAVID; http://david.abcc.ncifcrf.gov/) (32) for GO BP terms and KEGG pathway analyses. P<0.05 and gene counts >2 were considered to indicate a statistically significant difference.

PPI network construction

The Search Tool for the Retrieval of Interacting Genes database (STRING; http://www.bork.embl-heidelberg.de/STRING/) (33) provides comprehensive information on the functional interactions between DEGs and other genes by calculating their combined score. PPI pairs with a combined score of >0.7 as a cutoff value were identified. PPI networks were then constructed using Cytoscape software (version 2.6.3) (available at http://cytoscapeweb.cytoscape.org/) (34) based on the PPI pairs.

Statistical analysis

The correlation between the expression levels and the methylation levels of genes was analyzed by Pearson correlation coefficient test (35) using SPSS 19.0 software (IBM Corp., Armonk, NY, USA). The correlation coefficient takes a value between −1 and 1, where the value closer to −1 (or 1) indicates a negative (or positive) high association and closer to 0 indicates no association between the two variables.

Results

A total of 663 DEGs with log FC>1 and adjusted P<0.05 were screened using the Limma package of R, including 227 upregulated and 436 downregulated DEGs. The resultant heat map of DEGs is depicted in Fig. 1.

Figure 1.

Heat map of DEGs. Samples with the serial numbers GSM879175-GSM879180 represent normal controls, and the others represent OS cell lines. The left vertical strips represent DEGs in OS cell lines compared with normal controls, with yellow strips representing downregulated genes and blue strips representing upregulated ones. DEGs, differentially expressed genes; OS, osteosarcoma.

A total of 2,368 DMRs were identified using methyAnalysis with minimum P-values adjusted to <0.01, and 1,093 hyper- and 162 hypo-methylated regions were located in the gene promoter regions. Notably, there were 47 downregulated genes in the 1,093 hyper-methylated regions, while no genes were upregulated in the 162 hypo-methylated regions. Pearson correlation coefficient tests (35) was then used to calculate the correlation between the gene expression levels and the methylation levels of 47 downregulated genes. The correlation results demonstrated that the median, minimum and maximum values of the Pearson correlation coefficient were −0.4987, −0.7902 and −0.1573, respectively. The genes with the strongest correlations were intercellular adhesion molecule 3 (r=−0.7902), cytochrome c oxidase subunit VIIa polypeptide 1 (r=−0.7743), aminomethyltransferase (r=−0.7554), retinoic acid receptor responder (tazarotene induced) 1 (r=−0.7517) and Fc fragment of IgE, high affinity I, receptor for; gamma polypeptide (FCER1G; r=−0.7366).

Based on the microarray data for SNP profiling, 2,838 genes exhibited copy number gains and 415 exhibited copy number losses. Among them, 17 genes with copy number gains were upregulated, while 5 genes with copy number losses were downregulated. The genes with upregulated expression and copy number gains were as follows: ATPase family, AAA domain containing 2 (ATAD2), small nuclear ribonucleoprotein polypeptide A', ribosomal protein L7 (RPL7), cyclin-dependent kinase 4 (CDK4) and tubulin, α 1a (TUBA1A). The genes with downregulated expression and copy number losses were Rho GTPase activating protein 9 (ARHGAP9), chromosome 1 open reading frame 54, leiomodin 3 (fetal), major histocompatibility complex, class II, DO α (HLA-DOA) and Rh blood group, D antigen (RHD).

DAVID was used to perform functional enrichment analysis of the DEGs. The over-represented GO BP terms and significantly enriched KEGG pathways are summarized in Tables I and II, respectively. According to the enrichment results, the upregulated genes were mainly associated with carboxylic acid biosynthetic, organic acid biosynthetic and L-serine metabolic processes (Table I). In addition, the upregulated genes were significantly enriched in 6 KEGG pathways, including glycine, serine and threonine metabolism, ribosome and one carbon pool by folate (Table II). Downregulated genes were mainly associated with functions including the immune response, defense response and gas transport (Table I), and were significantly enriched in 20 KEGG pathways, including asthma, viral myocarditis and intestinal immune network for immunoglobulin (Ig) A production (Table II).

Table I.

The top 10 GO terms enriched among differentially expressed genes.

GO ID	Description	Counts	P-value
Upregulated
GO:0046394	Carboxylic acid biosynthetic process	13	1.33×10−8
GO:0016053	Organic acid biosynthetic process	13	1.33×10−8
GO:0006563	L-serine metabolic process	5	4.15×10−7
GO:0008652	Cellular amino acid biosynthetic process	7	5.86×10−6
GO:0006412	Translation	13	4.05 ×10−5
GO:0009069	Serine family amino acid metabolic process	5	7.82×10−5
GO:0009309	Amine biosynthetic process	7	8.49×10−5
GO:0009070	Serine family amino acid biosynthetic process	4	1.09×10−4
GO:0006564	L-serine biosynthetic process	3	2.37×10−4
GO:0006633	Fatty acid biosynthetic process	6	7.00×10−4
Downregulated
GO:0006955	Immune response	68	9.91×10−24
GO:0006952	Defense response	52	3.80×10−15
GO:0015669	Gas transport	9	3.36×10−9
GO:0009611	Response to wounding	38	4.05×10−9
GO:0032101	Regulation of response to external stimulus	20	7.43×10−9
GO:0048584	Positive regulation of response to stimulus	23	4.88×10−8
GO:0002504	Antigen processing and presentation of peptide or polysaccharide antigen via MHC class II	10	5.12×10−8
GO:0002684	Positive regulation of immune system process	23	5.67×10−8
GO:0050727	Regulation of inflammatory response	13	2.14×10−7
GO:0006954	Inflammatory response	26	2.54×10−7

GO ID represents the identification number of the enriched GO term. Description represents the name of the GO term. Counts represent the number of upregulated or downregulated target genes. GO, gene ontology.

Table II.

The significant KEGG pathways enriched by differentially expressed genes.

Description	Counts	%	P-value
Upregulated
hsa00260: Glycine, serine and threonine metabolism	5	3.125	6.26×10−04
hsa03010: Ribosome	6	3.75	0.005144
hsa00670: One carbon pool by folate	3	1.875	0.017709
hsa00100: Steroid biosynthesis	3	1.875	0.019904
hsa01040: Biosynthesis of unsaturated fatty acids	3	1.875	0.032434
hsa00450: Selenoamino acid metabolism	3	1.875	0.04415
Downregulated
hsa05310: Asthma	11	2.820513	9.14×10−09
hsa05416: Viral myocarditis	13	3.333333	1.71×10−06
hsa04672: Intestinal immune network for IgA production	11	2.820513	2.17×10−06
hsa05330: Allograft rejection	9	2.307692	1.21×10−05
hsa04514: Cell adhesion molecules (CAMs)	16	4.102564	1.45×10−05
hsa05332: Graft-versus-host disease	9	2.307692	2.26×10−05
hsa04940: Type I diabetes mellitus	9	2.307692	4.00×10−05
hsa05322: Systemic lupus erythematosus	13	3.333333	5.72×10−05
hsa04640: Hematopoietic cell lineage	12	3.076923	7.14×10−05
hsa05320: Autoimmune thyroid disease	9	2.307692	1.69×10−04

Description represents the identification number and the name of the KEGG pathway. Counts represent the number of upregulated or downregulated target genes. KEGG, Kyoto Encyclopedia of Genes and Genomes.

Notably, the enrichment results demonstrated that 47 downregulated genes in the hyper-methylated regions were significantly enriched in the ECM-receptor interaction pathway, and were associated with the regulation of cytokine production, regulation of tumor necrosis factor production and cell adhesion.

PPI network analysis

Based on the information from the STRING database, PPI networks of upregulated and downregulated genes were constructed. MAD2 mitotic arrest deficient-like 1 (degree=9), BUB1 mitotic checkpoint serine/threonine kinase (degree=8) and PDZ binding kinase (degree=7) were selected as the hub nodes in the PPI network of upregulated genes (Fig. 2). Notably, the upregulated genes, including ATAD2, RPL7 and CDK4, exhibited copy number gains and may interact with other upregulated genes (Fig. 2).

Figure 2.

Protein-protein interaction network of upregulated genes. Red nodes represent upregulated genes. Blue nodes represent upregulated genes with copy number gain.

In the PPI network of downregulated genes, TYRO protein tyrosine kinase binding protein (degree=37), immunoglobulin superfamily, member 6 (degree=20) and lysosomal protein transmembrane 5 (degree=19) with the highest degrees were hub nodes (Fig. 3). In addition, the hyper-methylated hub nodes were FCER1 G, leptin (LEP), and feline Gardner-Rasheed sarcoma viral oncogene homolog (FGR). Concurrently, ARHGAP9, HLA-DOA and RHD exhibited copy number losses and may interact with other downregulated genes.

Figure 3.

Protein-protein interaction network of downregulated genes. Green nodes represent downregulated genes. (A) Blue nodes represent downregulated genes with copy number loss. (B) Yellow nodes represent hyper-methylated downregulated genes.

Discussion

The cell line panel provides a valuable model system for analysis of genetic and epigenetic aberrations in OS (36). In the present study, a comprehensive bioinformatics approach was utilized for analysis of the effects of genome-wide changes in gene expression, DNA methylation and DNA copy number alterations in OS cell lines compared with the normal controls. The results demonstrated that a total of 47 downregulated genes were located in hyper-methylated regions, while no genes were upregulated in hypo-methylated regions, including FCER1G, FGR and LEP. In addition, downregulated genes including ARHGAP9 and HLA-DOA, exhibited copy number losses, and upregulated genes, including ATAD2 and CDK4, exhibited copy number gains. All these genes may act a key function in OS progression and merit additional discussion.

FCER1G, LEP and FGR were hyper-methylated in the present study. FCER1G, the γ subunit gene of the high-affinity receptor for IgE (FcεRI), is considered to serve a crucial function in upregulating FcεRI on atopic antigen-presenting cells (37). Demethylation of the FCER1G promoter results in the overexpression of FcεRI on monocytes of patients with atopic dermatitis (38). Epigenetic deregulation may serve an important function in cancer development (39). LEP encodes the leptin protein which is expressed in and secreted from mature primary cultures of human osteoblasts, which is of importance for osteoblastic cell growth and bone mineralization (40). LEP may increase bone mass by stimulating OS cell proliferation via activation of the phosphoinositide 3-kinase and mitogen-activated protein kinase (MAPK) signaling pathways (41). The FGR gene is a member of the Src family of protein tyrosine kinases. This gene functions as a negative regulator of cell migration and adhesion, and is triggered by the β-2 integrin signal transduction pathway (42). The chemokine receptor 4/stromal cell-derived factor 1 serves an important function in OS tumor progression via the regulation of cell migration and adhesion (43), suggesting that FGR may be involved in OS development via regulating cell migration and adhesion. In the present study, these hyper-methylated genes were revealed to be downregulated, thus, it was hypothesized that these genes may be involved in OS progression via decreased expression following hyper-methylation. The results may assist in understanding the epigenetic regulation of specific genes in OS tumor development.

In addition, OS is a copy number alteration-rich malignant bone tumor, and candidate genes with copy number changes are being identified in OS (44). In the present study, DEGs including ARHGAP9 and HLA-DOA, exhibited copy number losses while ATAD2 and CDK4 exhibited copy number gains. ARHGAP9, which is a novel MAP kinase docking protein, interacts with mitogen-activated protein kinase 1 (ERK2) and p38α (45). Silencing ERK2 in human U2OS cells may inhibit the expression and function of glycoprotein 130, which serves a pivotal function in cancer and inflammation (46). In addition, HLA-DOA is a key molecule in the antigen processing and presentation pathway, and this pathway has also been suggested to be involved in OS progression via downregulated expression of HLA-DOA (47,48). In addition, ATAD2 is highly expressed and genetically amplified in several types of human cancer (49). ATAD2 binds to the v-Myc avian myelocytomatosis viral oncogene homolog (c-Myc) oncogene and stimulates its transcriptional activity (49). A previous study demonstrated that overexpression of c-Myc may promote OS cell invasion via the activation of the mitogen activated protein kinase kinase-extracellular signal-regulated kinase pathway (50). ATAD2 is also confirmed to exhibit prognostic significance in high-grade OS (51). The amplification and overexpression of CDK4 tends to be associated with improved prognosis in low-grade OS (52). 12q13-14 CDK4 amplicons are frequently observed in OS (53). CDK4 and other CDK inhibitors are regarded as promising anticancer agents in cancer treatment (54). Therefore, these DEGs may serve important functions in the development and progression of OS. The results of the present study are consistent with these data, suggesting that copy number alterations of key genes may be associated with OS.

In conclusion, hyper-methylation of FCER1G, LEP and FGR is observed in OS, suggesting that epigenetic alterations of these specific genes may act crucial functions in OS development. In addition, copy number alterations of these DEGs, including ARHGAP9, HLA-DOA, ATAD2 and CDK4, may also contribute to OS progression. These results indicate that genetic and epigenetic alterations are important mechanisms involved in OS, and these DEGs may serve as candidate targets for the diagnosis and treatment of this disease. However, no experimental validation and the relatively small sample size are the limitations of the present study. The results require additional validation.