Identification and analysis of key genes in osteosarcoma using bioinformatics.

Osteosarcoma (OS) is an invasive malignant neoplasm of the bones. The present study identified and analyzed key genes associated with OS. Expression profiling of the dataset GSE49003, which included 6 metastatic and 6 non-metastatic OS cell lines and was obtained from the Gene Expression Omnibus, was performed. Following data preprocessing, the differentially expressed genes (DEGs) were selected using the limma package in R. Subsequently, bidirectional hierarchical clustering using the pheatmap package in R and an unpaired Students' t-test was performed for the DEGs. Based on the Search Tool for the Retrieval of Interacting Genes database and Cytoscape software, a protein-protein interaction (PPI) network for the DEGs was constructed. Using Database for Annotation, Visualization and Integrated Discovery software and the Kyoto Encyclopedia of Genes and Genomes Orthology Based Annotation System server, functional and pathway enrichment analyses were performed for the DEGs corresponding to the proteins of the network. In addition, the transcription factors (TFs) and CpG islands of the gene promoter were searched for using the TRANSFAC database and CpG Island Searcher software, respectively. A total of 323 DEGs were identified between the metastatic and non-metastatic samples. In the PPI network, upregulated epidermal growth factor receptor (EGFR) exhibits a high degree and was therefore highly interconnected with other proteins. Enrichment analysis revealed that EGFR was enriched in cytoskeleton organization, organic substance response and the signaling pathway of focal adhesion. The TFs early growth response 1, nuclear factor-κB complex subunits, peroxisome proliferator activated receptor α, signal transducer and activator of transcription 3 and MYC proto-oncogene were identified in the EGFR promoter region. Furthermore, multiple CpG islands, starting from the 400 bp of the EGFR promoter sequence, were predicted. Methylated modification of the CpG islands in the EGFR promoter may help to regulate EGFR expression. The TFs identified in the EGFR promoter may function in the progression of OS.

Introduction

Osteosarcoma (OS) is an invasive malignant neoplasm of the bones that originates from primitive transformed cells (1). OS is often derived from the metaphyseal region of tubular long bones, with 10, 19, and 42% of OS cases occurring in the humerus, tibia, and femur, respectively (2). OS is the most common type of bone cancer among children and teenagers (3). In pediatric patients, OS ranks eighth among all types of cancer incidence, accounting for ~20% of all primary bone cancers and 2.4% of all malignant tumors (4). Therefore, assessing the key genes associated with OS is crucial.

The downregulation of naked cuticle homolog 2 and alterations in its associated signaling pathways contribute to tumor growth and metastasis in patients with OS (5). Ribophorin II may represent a promising therapeutic target for OS, as silencing it promotes the regulation of lethal OS phenotypes (6). Cluster of differentiation 276 molecules belongs to the B7 family of costimulatory molecules and is associated with metastasis and survival in patients with OS, and its expression in OS cells may serve to control invasive malignancy and tumor immunity (7). By inhibiting epithelial cadherin expression and inducing the epithelial-mesenchymal transition (EMT) in OS cells, the overexpression of the transcription factor (TF) snail family transcriptional repressor 1 (SNAI1) serves a function in tumor progression by promoting invasion and metastasis (8). Overexpression of EMT-associated TFs, including SNAI1, twist family bHLH transcription factor 1 and zinc finger E-box binding homeobox 1 is involved in the progression of osteosarcoma, and molecules or chemotherapeutic agents targeting these TFs may be used to treat OS by helping to control metastasis (9). Furthermore, Al-Romaih et al (10) demonstrated that a significant loss of 5′ CpG island methylation in growth arrest and DNA damage inducible α is associated with the induction of apoptosis in OS cells. However, the mechanisms underlying OS cell apoptosis remain to be fully understood.

In 2015, Endo-Munoz et al (11) identified molecules that may promote OS metastasis through multiomics analysis, confirmed this function in vitro and in vivo, and demonstrated that the plasminogen activator, urokinase (UPA)/UPA receptor axis contributed to OS metastasis. Since different procedures may yield different results, the present study further analyzed the data from the Endo-Munoz et al study (11) using bioinformatics. Initially, the differentially expressed genes (DEGs) in the metastatic OS cell lines were identified. Secondly, bidirectional hierarchical clustering and an unpaired Student's t-test were performed to confirm the DEGs between the metastatic and the non-metastatic samples. Thirdly, a protein-protein interaction (PPI) network for the DEGs was constructed and enrichment analysis was conducted for the DEGs corresponding to the proteins of this network, to identify key genes. Finally, the TFs and CpG islands of the promoter regions of key genes were identified.

Materials and methods

Data source

The microarray GSE49003 dataset (11) was obtained from the database of the Gene Expression Omnibus (www.ncbi.nlm.nih.gov/geo) (12), and was sequenced on the GPL6947 Illumina HumanHT-12 V3.0 expression beadchip (Illumina, Inc., San Diego, CA, USA). GSE49003 included 3 metastatic OS KRIB cell lines, 3 metastatic OS KHOS cell lines, 3 non-metastatic OS U2OS cell lines, and 3 non-metastatic OS HOS cell lines. These OS cell lines were provided by the Basil Hetzel Institute of the University of Adelaide (Adelaide, Australia). The cells were cultured in Dulbecco's modified Eagle's medium supplemented with 10% fetal calf serum and pregnancy-specific glycoproteins in a 5% CO2 atmosphere for 24 h at 37°C until 80% confluent (11).

Data preprocessing and differential expression analysis

Based on the platform annotation information of GSE49003, probes were transformed into their corresponding gene symbols. For probes that corresponded to the same gene, their values were averaged to obtain the gene expression value. Subsequently, log2 transformation and normalization were performed for the data. The DEGs between non-metastatic and metastatic samples were identified using the limma package in R (www.bioconductor.org/packages/release/bioc/html/limma.html) (13). Using the Benjamini-Hochberg method (14) via the multtest package in R, the P-values were adjusted according to the false discovery rate (FDR; thresholds, FDR <0.05 and |log2 fold change| >1). Following the extraction of the gene expression values in each sample, bidirectional hierarchical clustering (15) based on Euclidean distance (16) was performed using the pheatmap package in R (https://cran.r-project.org/web/packages/pheatmap/index.html) (17). To determine whether the DEGs distinguished metastatic and non-metastatic samples, an unpaired Student's t-test (18) was used to compare the expression of the DEGs between the samples. P<0.05 was considered to indicate a statistically significant difference.

PPI network analysis

The Search Tool for the Retrieval of Interacting Genes (STRING; http://string-db.org/) includes verified and predicted interactions (19). Cytoscape software (version 3.2.0; Institute for Systems Biology, Seattle, WA, USA) was used to construct a unified framework by integrating gene expression data with biomolecular interaction networks (20). The interactions among the DEGs were searched for using the STRING database (version 10.0) (19), with a combined score >0.8 as the threshold. Subsequently, the PPI network was visualized using Cytoscape software (version 3.2.0) (20). The number of interactions associated with a given node corresponded to the degree for the node. In addition, the nodes that exhibited increased degrees were hub nodes.

Enrichment analysis for the DEGs corresponding to the proteins of the PPI network

Gene Ontology (GO; http://www.geneontology.org/) includes the molecular functions, biological processes and cellular components associated with genes and their products (21). The Kyoto Encyclopedia of Genes and Genomes (KEGG; www.genome.jp/kegg) integrates chemical, genomic and systemic functional information on genes, and may be applied in pathway mapping (22). Using Database for Annotation, Visualization and Integrated Discovery software (version 6.7; david.abcc.ncifcrf.gov) (23), the DEGs corresponding to the proteins of the PPI network underwent functional enrichment analysis, with an FDR threshold of <0.05. Furthermore, the KEGG Orthology-Based Annotation System server (kobas.cbi.pku.edu.cn/home.do) (24) was used to perform pathway annotation and enrichment analysis for the DEGs corresponding to the proteins of the PPI network, with P<0.05 was considered to indicate a statistically significant difference.

Search for the TFs and CpG islands of the promoter regions of key genes

TRANSFAC (http://www.gene-regulation.com/pub/databases.html) is a database that integrates eukaryotic transcriptional regulatory information, including information on TFs, transcription regulation relationships and TF binding sites (25). Using the TRANSFAC database (version 11.2) (25), the TFs in the promoter regions of key genes were searched for. CpG islands are present in the majority of mammalian gene promoters and are characterized by an increased frequency of non-methylated CpG dinucleotides (26). Abnormal CpG island hypermethylation is common in certain types of cancer and may promote tumor progression by inhibiting tumor suppressor gene expression (27). The promoter sequences of the key genes were downloaded from Ensembl (http://pre.ensembl.org) (28), and their CpG islands were predicted using CpG Island Searcher software (cpgislands.usc.edu) (29). The thresholds were as follows: GC >50%, CpG observed/expected value >0.6, and CpG island length >200 bp.

Results

DEG analysis

A total of 323 DEGs were identified between the metastatic and non-metastatic samples. Of these, 189 genes were upregulated and 134 were downregulated. Bidirectional hierarchical clustering indicated that the DEGs differentiated the metastatic samples from the non-metastatic samples (Fig. 1A). The t-test revealed that the DEGs were significantly differentially expressed and more likely to be upregulated than downregulated (Fig. 1B).

Figure 1.

DEGs between the metastatic and non-metastatic samples. (A) Hierarchical clustering tree. Orange indicates the upregulated DEGs and blue indicates the downregulated DEGs. (B) Box plot for the DEGs. DEG, differentially expressed gene.

PPI network construction and enrichment analysis

The PPI network consisted of 80 nodes (42 corresponding to the products of upregulated genes; 38 corresponding to the products of downregulated genes) and 93 interactions (Fig. 2). Upregulated expression of epidermal growth factor receptor (EGFR, degree=12) was associated with an increased degree. For the DEGs corresponding to the proteins of the PPI network, a total of 7 GO terms were enriched, including cytoskeleton organization (FDR=0.003084; associated with EGFR), response to endogenous stimulus (FDR=0.008481) and response to organic substance (FDR=0.009596; associated with EGFR; Table I). Only the signaling pathways of focal adhesion (P=0.006078; associated with EGFR) and extracellular matrix-receptor interaction (P=0.013264) were enriched for the DEGs corresponding to the proteins of the PPI network (Table II).

Figure 2.

PPI network constructed for the differentially expressed genes between the metastatic and non-metastatic samples. Red and blue nodes denote upregulated and downregulated genes, respectively. PPI, protein-protein interaction.

Table I.

Enriched functions for the differentially expressed genes, corresponding to the proteins of the protein-protein interaction network.

Term	Count	P-value	False discovery rate	Gene symbol
GO:0007010 - Cytoskeleton organization	10	3.29×10−5	0.003084	EGFR, BMP4, MSX1, COL13A1, HOXA5, CTGF, JUND, TGFBR2, COL12A1, COL1A1
GO:0009719 - Response to endogenous stimulus	11	3.63×10−5	0.008481	BMP4, DUSP1, CCND2, PPARG, TGFBR2, HMGCS1, RHOQ, CA2, COL1A1, FOSL1, CCNA2
GO:0010033 - Response to organic substance	16	2.64×10−6	0.009596	EGFR, BMP4, TGFBR2, PPARG, HMGCS1, NFKBIA, RHOQ, MSX1, CD44, DUSP1, CCND2, ID1, CA2, COL1A1, FOSL1, CCNA2
GO:0009725 - Response to hormone stimulus	10	9.72×10−5	0.01414	BMP4, DUSP1, CCND2, PPARG, TGFBR2, RHOQ, CA2, COL1A1, FOSL1, CCNA2
GO:0042127 - Regulation of cell proliferation	14	1.47×10−4	0.015592	BMP4, EGFR, S100A6, NDN, TGFBR2, PPARG, NFKBIA, LAMA1, MSX1, RAC2, CCND2, FOSL1, CCNA2, PLAU
GO:0006928 - Cell motion	11	1.37×10−4	0.015959	ACTB, LAMA1, NDN, CD44, ID1, CTGF, CALD1, TPM1, PLAU, PLAUR, ARHGDIB
GO:0006357 - Regulation of transcription from RNA polymerase II promoter	13	2.77×10−4	0.022896	BMP4, KLF6, FOXA2, PPARG, RHOQ, NFKBIA, MSX1, HEY1, ID1, JUND, IRF8, FOSL1, NKX2-2

Table II.

Pathways enriched for the differentially expressed genes corresponding to the proteins of the protein-protein interaction network.

ID	Signaling pathway	Count	P-value	Gene symbol
hsa04510	Focal adhesion	8	0.006078	EGFR, ACTB, LAMA1, RAC2, CCND2, TNC, COL6A3, COL1A1
hsa04512	Extracellular matrix-receptor interaction	5	0.013264	LAMA1, CD44, TNC, COL6A3, COL1A1

Using the TRANSFAC database, the TFs in the promoter regions of EGFR were searched. Finally, the TFs early growth response 1 (EGR1), nuclear factor (NF)-κB complex subunits, peroxisome proliferator activated receptor α (PPARA), signal transducer and activator of transcription 3 (STAT3) and MYC proto-oncogene (MYC) were identified in the promoter region of EGFR (Fig. 3). The promoter sequence of EGFR was downloaded, and multiple CpG islands, starting from the 400 bp of the promoter sequence, were predicted (Fig. 4).

Figure 3.

Transcription factors in the promoter region of epidermal growth factor receptor. Egr, early growth response; NF, nuclear factor; PPAR, peroxisome proliferator activated receptor; STAT, signal transducer and activator of transcription; C-Myc, MYC proto-oncogene.

Figure 4.

Prediction map of the CpG islands in the promoter sequence of epidermal growth factor receptor. (A) Obs/exp CpG ratio, the red line presents Obs/exp CpG=0.6. (B) G and C content, the red line presents G and C content=60%. The threshold was set as follows: CpG island length >200 bp, G and C content >50%, Obs/exp CpG >0.6. Obs, observed; exp, expected.

Discussion

In the present study, a total of 323 DEGs were identified between the metastatic and non-metastatic samples. Of these, 189 genes were upregulated and 134 were downregulated. Bidirectional hierarchical clustering and a t-test indicated that the DEGs distinguished the metastatic and non-metastatic samples, and were significantly differentially expressed. In the PPI network constructed for the DEGs, upregulated EGFR exhibited a higher degree and therefore was identified to be highly interconnected with other proteins. Enrichment analysis revealed that EGFR was enriched in cytoskeleton organization, organic substance response, and the signaling pathway of focal adhesion. Expression and genomic gains of EGFR, and deletions of the phosphatase and tensin homolog, are common in OS (30). In addition, the EGFR signaling pathway may be inhibited by microRNA (miR)-143 by regulating matrix metallopeptidase 9 (MMP9) expression in OS; therefore, EGFR, MMP9 and miR143 may serve as therapeutic targets for inhibiting OS invasion (31). The nuclear expression of erb-b2 receptor tyrosine kinase 4 (ERBB4)and the cell surface expression of ERBB2 and EGFR may contribute to OS progression, indicating that therapy targeting the expression of these genes may benefit patients with OS (32). EGFR expression has been previously suggested to serve a critical function in the progression of OS, is increased in OS cell lines, and its downstream signaling molecule MAPK1 serves as a prognostic factor in patients with OS (33). Mammalian CpG island methylation may result in transcriptional silencing, and 3′ CpG island methylation serves a function in developmental gene regulation (34). Multiple CpG islands starting from the 400 bp of the EGFR promoter sequence were predicted, and these regions possessed an increased likelihood of undergoing methylation, compared with regions <400 bp, which may serve a function in regulating EGFR expression.

Using the TRANSFAC database, the TFs EGR1, NF-κB complex subunits, PPARA, STAT3 and MYC were identified in the promoter region of EGFR. EGR1 may be expressed in normal bone cells following stimulation with neurotransmitters, hormones and growth factors, which serve functions in bone formation (35). Lysophosphatidic acid-induced upregulation of EGR1 expression is accompanied by a time-dependent decrease in periostin expression in MG-63 osteosarcoma cells, and therefore EGR1 may influence periostin expression in these cells (36). Downregulating solute carrier family 16 member 1 (SLC16A1) expression has been demonstrated to regulate the NF-κB complex subunits signaling pathway, which may result in tumor growth and metastasis inhibition in OS, while increased SLC16A1 expression has been revealed to predict poor overall survival in patients with OS (37). Since PPAR agonists may induce apoptosis and inhibit proliferation in OS cells, effective PPAR agonists may be used as adjuvant therapeutic drugs for patients with OS (38). The TF STAT3 has been reported to serve functions in tumor cell survival, growth, differentiation, metastasis, apoptosis and drug resistance (39,40). In multidrug resistant OS cells, the STAT3 signaling pathway is activated, and STAT3 phosphorylation and nuclear translocation are inhibited by the synthetic oleanane triterpenoid CDDO-Me, which may induce apoptosis (41). Dihydrofolate reductase, MYC and other prognostic markers may be used to identify patients with early stage OS that exhibit an increased risk of adverse outcomes (42). Furthermore, the protein expression of MYC, p53, BCL2, and the apoptotic index may be used to predict the progression and prognosis of OS, and for selecting an optimal treatment (43). Therefore, the TFs EGR1, NF-κB complex subunits, PPARA, STAT3 and MYC identified in the promoter region of EGFR may also function in the progression of OS.

To conclude, the present study identified a total of 323 DEGs between the metastatic and non-metastatic samples. The methylation of the CpG islands, starting from the 400 bp of the EGFR promoter sequence, may be used to regulate EGFR expression. Furthermore, the TFs EGR1, NF-κB complex subunits, PPARA, STAT3 and MYC identified in the promoter region of EGFR may be associated with the progression of OS. However, the results of the present study, predicted by bioinformatics, require further study for validation.