This article is mentioned in:

Introduction

Osteosarcoma (OS), also termed bone sarcoma, originates from bone and particularly from the mesenchymal stem cell lineage (1). OS, the most common bone tumor, is highly aggressive and usually has poor prognosis (2). Additionally, OS primarily affects adolescents and children and ~60% of neoplasms occur in patients under the age of 20 (3,4). Current treatment frequently involves a combination of surgery and chemotherapy; however, OS still leads to a high mortality and morbidity, particularly in children and adolescents (1).

Currently, considerable progress has been made in identifying the critical factors in the development and progression of OS, including genes, pathways and copy number variants (CNVs) (5). Alterations of tumor suppressor gene expression including protein kinase, cAMP-dependent, regulatory, and type I a and deregulation of major signaling pathways such as the wingless-type MMTV integration site family, transforming growth factor-β, Notch and sonic hedgehog have been previously associated with OS (6,7). It has also been previously demonstrated that OS development is dependent on loss of P53 and enhanced by loss of retinoblastoma 1 (RB1) (8). CNVs are DNA segments 1 kb in length which are present in a variable population frequency in the genome (9). During the 1990s, CNVs with duplications and deletions were expressed as an inducement of a quantity of single gene disorders (10). Various differentially expressed genes (DEGs) and several candidate CNVs in OS have been identified to be involved in the development of OS by analyzing the microarray data and high-resolution single nucleotide polymorphism (SNP)/CNV arrays (11,12). However, frequently only one of these approaches has been used in previous studies to identify the candidate molecule, and the molecular mechanism of OS remains to be elucidated (9–11).

The proto-oncogene MET protein, a hepatocyte growth factor receptor, encodes tyrosine-kinase activity (13), which has been revealed to be aberrantly expressed in OS and closely associated with cancer (14–16). Therefore, overexpression of the MET oncogene may convert human primary osteoblasts (HOB) into OS cells.

The present study extracted the transcriptional and CNV profiles from Gene Expression Omnibus (GEO) database. The differentially expressed genes (DEGs) and CNVs were screened in the transcriptional profile of MET-HOB cells, which were previously turned into OS cells by lentiviral vector (LV)-driven overexpression of the MET oncogene. Subsequently, the shared genes in the expression and CNV profiles were analyzed. The present study obtained a series of candidate markers in OS and may provide the foundation for treatment of OS.

Materials and methods

Microarray and CNV data

Transcriptional profile (ID: GSE28256) was downloaded from the GEO database (www.ncbi.nlm.nih.gov/geo/) which was based on the platform of GPL6098 (Illumina humanRef-8 v1.0 expression beadchip) (17). The dataset contained 15 samples, including 6 HOB cell lines and 9 MET-HOBs clones, which were previously turned into osteosarcoma cells by over-expression of MET oncogene driven by a LV. The CNV data were extracted from the GSE32964 dataset in the GEO database, which included 36 samples for detecting SNP and 32 samples for CNV. A total of 32 CNV samples of OS tumor tissues based on the platform of GPL6985 (Illumina HumanCNV370-QuadV3 DNA Analysis BeadChip) were analyzed in the present study.

Data preprocessing

The probe-level data of the transcriptional profile were initially converted into expression values. Probes that mapped with the gene names labeled in the annotation platform were transformed using log2 and normalized using preprocess-Core package in R version 2.9.0. According to the annotation platform, the values of probes corresponding to the same transcript were averaged and then defined as the final expression value of a transcript. The PennCNV tool (version 2014 May 07; http://penncnv.openbioinformatics.org) was used in the subsequent processing of data, the profile of CNV samples was converted into specific format for PennCNV, which contained log R Ratio: LRR and B Allele frequency: BAF. In addition, to investigate the differences among samples, a heatmap was generated to compare their expression values using the Gplots package in R.

Identification of DEGs and CNV

A Student's t-test was conducted on the gene expression values between testing and control samples. The Linear Models for Microarray Data (LIMMA) package was used to normalize the data and identify the DEGs in MET-HOB samples compared with control samples using cut-offs of P<0.05 and |log2 fold-change (FC)|>2. Additionally, detect_cnv.pl in PennCNV was applied to select CNV areas and CNVs were identified with the cut-off criteria of >30% overlap within the cases. Genes shared in CNVs and DEGs were identified as critical genes associated with the development of OS.

Functional enrichment of DEGs

The Database for Annotation, Visualization and Integrated Discovery (DAVID) provides numerous comprehensive functional annotation which contributes to the understanding of the biological meanings behind abundant genes (18). Gene Ontology (GO) and the Kyoto Encyclopedia of Genes and Genomes (KEGG) analyses have become commonly used approaches for functional and pathway studies of large-scale genomic or transcription data, respectively (19). Therefore, they were used in the present study. Next, DAVID was used to screen the enriched GO terms and KEGG pathways in the DEGs. P<0.05 was used as a cut-off criterion.

Results

Data preprocessing and DEGs screening

In the present study, 24,350 probes were detected in the original data and 21,454 non-redundant genes were obtained following data preprocessing. The raw data in all samples have been normalized (Fig. 1). A total of 1601 DEGs were screened out in MET-HOBs compared with controls. Among these genes, 784 were upregulated in MET-HOBs and 817 were downregulated. The hierarchical clustering analysis revealed a clearly distinct expression of all DEGs between MET-HOBs and HOBs (Fig. 2).

Figure 1. Boxplots for genes in all samples (A) before and (B) after normalization. Red indicates MET proto-oncogene-transformed HOB samples and blue boxes indicate HOB samples. HOB, human primary osteoblasts.

Figure 2. Hierarchical cluster analysis for the differentially expressed genes in MET-HOB and HOB samples. Cases 1–9 indicate MET-HOB samples and control 1–6 represent HOB samples. The gradual color change from green to red represents the changing process from upregulation to downregulation. HOB, human primary osteoblasts; MET-HOB, MET proto-oncogene-transformed human primary osteoblasts.

Function enrichment of DEGs

In order to identify the functions of the DEGs, they were performed GO (P<0.01) and KEGG (P<0.01) enrichment analyses. The results indicated that 344 GO terms were obtained and the top 10% terms were listed in Table I, such as extracellular region (P=2.68E-24), extracellular matrix (ECM; P=4.08E-24) and proteinaceous extracellular matrix (P=4.79E-24). Besides, 14 KEGG pathways were obtained and most of them were related to cancers, such as hsa05200: pathways in cancer (P=4.12E-09), hsa04512: ECM-receptor interaction (P=2.41E-08) and hsa05222: small cell lung cancer (P=1.06E-05; Table II).

Table I. The top 10% enriched GO terms for DEGs.

Table I.

The top 10% enriched GO terms for DEGs.

Category	GO ID	GO name	Gene number	P-value
CC	GO:0044421	Extracellular region part	134	2.68×10−24
CC	GO:0031012	Extracellular matrix	74	4.08×10−24
CC	GO:0005578	Proteinaceous extracellular matrix	71	4.79×10−24
BP	GO:0001501	Skeletal system development	60	4.57×10−15
MF	GO:0019838	Growth factor binding	32	2.32×10−14
BP	GO:0001944	Vasculature development	50	1.22×10−13
BP	GO:0001568	Blood vessel development	49	1.88×10−13
MF	GO:0005201	Extracellular matrix structural constituent	28	2.11×10−13
BP	GO:0042127	Regulation of cell proliferation	98	4.21×10−12
BP	GO:0006260	DNA replication	40	8.03×10−12
BP	GO:0007049	Cell cycle	96	1.09×10−11
BP	GO:0051726	Regulation of cell cycle	55	1.23×10−11
BP	GO:0022403	Cell cycle phase	63	1.57×10−11
CC	GO:0005576	Extracellular region	179	1.87×10−11
CC	GO:0044427	Chromosomal part	57	2.60×10−11
BP	GO:0051270	Regulation of cell motion	39	5.52×10−11
BP	GO:0006259	DNA metabolic process	70	8.61×10−11
CC	GO:0044420	Extracellular matrix part	28	1.55×10−10
BP	GO:0040012	Regulation of locomotion	38	1.89×10−10
BP	GO:0051301	Cell division	49	1.93×10−10
BP	GO:0030334	Regulation of cell migration	35	3.13×10−10
BP	GO:0007155	Rell adhesion	85	5.11×10−10
BP	GO:0022610	Biological adhesion	85	5.66×10−10
CC	GO:0005615	Extracellular space	79	6.44×10−10
BP	GO:0006928	Cell motion	65	6.97×10−10
BP	GO:0022402	Cell cycle process	73	7.17×10−10
CC	GO:0005694	Chromosome	60	1.09×10−9
BP	GO:0000278	Mitotic cell cycle	54	2.53×10−9
BP	GO:0000279	M phase	49	7.98×10−9
BP	GO:0048514	Blood vessel morphogenesis	37	1.03×10−8
BP	GO:0065004	Protein-DNA complex assembly	23	1.29×10−8
CC	GO:0005581	Collagen	14	1.97×10−8
BP	GO:0030198	Extracellular matrix organization	24	3.69×10−8
BP	GO:0008283	Cell proliferation	57	4.64×10−8

[i] GO, Gene Ontology; DEGs, differentially expressed genes; CC, cellular component; MF, molecular function; BP, biological process.

Table II. Enriched Kyoto Encyclopedia of Genes and Genomes pathways for DEGs.

Table II.

Enriched Kyoto Encyclopedia of Genes and Genomes pathways for DEGs.

KEGG ID	Pathway name	Gene number	P-value
hsa05200	Pathways in cancer	53	4.12×10−9
hsa04512	ECM-receptor interaction	23	2.41×10−8
hsa05222	Small cell lung cancer	19	1.06×10−5
hsa03030	DNA replication	12	1.72×10−5
hsa04510	Focal adhesion	30	8.78×10−5
hsa04115	p53 signaling pathway	15	1.62×10−4
hsa05210	Colorectal cancer	16	4.93×10−4
hsa05218	Melanoma	14	9.09×10−4
hsa05219	Bladder cancer	10	1.79×10−3
hsa00980	Metabolism of xenobiotics by cytochrome P450	12	2.20×10−3
hsa05215	Prostate cancer	15	2.68×10−3
hsa05217	Basal cell carcinoma	11	3.69×10−3
hsa05214	Glioma	11	9.89×10−3

[i] DEGs, differentially expressed genes.

Identification of CNVs

A total of 1,313 chromosome regions were identified and 678 genes (239 duplications and 439 deletions) were obtained, which were spread among 22 pairs of autosomes (Fig. 3). Then these CNVs were checked for overlap with the DEGs. Finally, 12 genes were identified in both in CNVs and DEGs, including the six upregulated genes cadherin 18 (CDH18), spectrin β, erythrocytic (SPTB), ciliary rootlet coiled-coil, rootletin pseudogene 2 (CROCCP2), β-1,4-N-acetyl-galactosaminyltransferase 1 (B4GALNT), G protein regulated inducer of neurite outgrowth 1 (GPRIN1) and growth factor independent 1 (GFI1). A total of six downregulated genes were identified, including laminin subunit α 1 (LAMA1), EH domain binding protein 1-like 1 (EHBP1L1), cathepsin Z (CTSZ), WNK lysine deficient protein kinase 1 (PRKWNK1), glutathione S-transferase µ 2 (GSTM2) and microsomal glutathione S-transferase 1 (MGST1) (Table III).

Figure 3. Genome distributions of CNVs in MET proto-oncogene-transformed human primary osteoblast samples. The PennCNV tool was used to select CNV areas. CNVs were identified with the cut-off criteria of >30% overlap within the dataset. Red represents deletion and purple represents duplication. CNVs, copy number variations.

Table III. CNV-driven genes.

Table III.

CNV-driven genes.

Gene	logFC
LAMA1	−2.583
EHBP1L1	−5.75702
SPTB	5.224308
CTSZ	−7.39063
CROCCP2	2.491733
B4GALNT	3.6299
CDH18	2.227242
GPRIN1	2.361558
PRKWNK1	−2.75581
GFI1	3.577933
GSTM2	−2.62289
MGST1	−2.30508

[i] CNV, copy number variations; FC, fold-change.

Discussion

OS is a universally fatal disease, due to the rapid growth, high local aggressiveness, and metastasizing potential (20). Numerous DEGs and regulatory relationships between transcription factors and DEGs in OS have been identified using microarray data (12). Additionally, susceptibility genes associated with OS were also reported by analyzing SNP/CNV arrays (11). However, the underlying molecular mechanism of OS remains to be elucidated. In the present study, 1,601 DEGs were identified, including 784 upregulated and 817 downregulated DEGs and CNVs in 678 genes (239 duplications and 439 deletions) were observed in MET-HOBs samples when compared with controls. By analyzing the transcriptional profile and SNP/CNV arrays, CDH18, LAMA1, SPTB, CROCCP2, B4GALNT, GPRIN1, GFI1, EHBP1L1, CTSZ, PRKWNK1, GSTM2 and MGST1 were identified as CNV-driven DEGs.

The DEGs obtained in the current study suggested that several genes such as E2F transcription factor 1 (E2F1) and 2 (E2F2), retinoblastoma 1 (RB1) and cyclin D1 (CCND1) were involved in various pathways. E2F1 and E2F2, members of the E2F family of transcription factors, were upregulated in the MET-HOBs samples. E2F proteins regulate the transcription of genes required for DNA synthesis (21). The E2F family has an important role in cell cycle regulation and action of tumor suppressor proteins, and is also a target of the transforming proteins of small DNA tumor viruses (22,23). Additionally, the RB protein has been previously identified to bind to E2F transcription factors (24). It is evident that the RB/E2F pathway is very important in regulating the initiation of DNA replication and the pathway is disrupted in the majority of human cancers (25). CCND1, is part of the highly conserved cyclin family, is a nuclear protein required for cell cycle progression in G1 phase (26). It has been previously reported that CCND1 has an important role in the regulation of OS cell proliferation (26). Consistently, the findings of the present study revealed that those genes were involved in several cell cycle-associated GO terms, such as cell cycle, cell cycle phase, regulation of cell cycle and cell division, and cancer-associated KEGG pathways, including pathways in cancer, small cell lung cancer and melanoma. Therefore, the present study is reliable and may suggest that the screened DEGs such as E2F1, E2F2, RB1 and CCND1 are closely associated with the cell cycle and cell division of OS.

CNVs such as deletions, duplications and amplifications across the whole genome may contribute to OS tumorigenesis (27). It is of note that CNVs in cyclin-dependent kinase inhibitor 2A (CDKN2A), sex determining region Y-box 6 (SOX6) and phosphatase and tensin homolog (PTEN) were associated with Ewing sarcoma (28). The trail of CDKN2A/B locus was detected in OS cell lines (29), whereas two SNPs in the SOX6 gene were identified to be associated with both hip bone mineral density and body mass index in Caucasians (30). In addition, copy number losses in PTEN were common events in OS (31). In the present study, further analysis identified 12 CNV-driven genes in MET-HOBs samples, such as CDH18 (upregulated) and LAMA1 (downregulated), which were associated with cell adhesion. Since the 1990s, many cadherins and cadherin-associated proteins had been identified and implicated in cancers as candidate tumor suppressors or proto-oncogenes (32). Deregulation of cadherin-catenin complexes may contribute to tumor development by influencing the adhesion of epithelial cells (33). CDH18 is a member of the cadherin superfamily that mediates calcium-dependent cell-cell adhesion (34). Although, no previous studies have not identified a direct association between CDH18 and OS, it has been revealed that an Exon 2 deletion of CDH18 may be associated with human colorectal cancer and CDH18 may act as novel candidate gene involved in colorectal cancer predisposition (35). In the current study, CDH18 was upregulated and CNVs were detected in CDH18 of MET-HOBs samples; therefore, CDH18 may have a role in cell-cell adhesion of OS. Additionally, LAMA1, also termed EHS laminin, was downregulated in MET-HOBs samples when compared with controls. Laminin is a complex glycoprotein and is considered to control the attachment, migration and organization of cells during embryonic development by interacting with other ECM components (36,37). Additionally, it has been previously reported that metadherin, as a laminin receptor, has an important role in controlling tumorigenesis and metastasis in many human cancers (38–40). Metadherin, a type II membrane protein in OS cells, may enhance cell invasion by regulating cell adhesion to the ECM through interaction with laminin (41). Therefore, LAMA1 may be involved in cell adhesion and cell-cell interactions in OS.

It is of note that there were some limitations in the present study. Only the OS-associated CNVs and DEGs were identified, whereas the transcription factors and protein-protein interaction network remain to be determined. The current findings were obtained by bioinformatics analysis and the corresponding validations were not performed. Therefore, future studies should involve in performing experiments such as reverse transcription-quantitative polymerase reaction and western blotting to validate the CNVs and DEGs identified.

In conclusion, a series of DEGs were identified to be associated with cell cycle and cell division of human OS, specifically E2F1, E2F2, RB1 and CCND1. Additionally, 12 CNV-driven DEGs were obtained and the cell adhesion-associated genes such as CDH18 and LAMA1 may contribute to OS cell-cell adhesion. These genes may act as alternative diagnosis and/or therapeutic markers for patients with OS. The present study developed the current understanding about the etiology of OS and provided the foundation for the development of novel treatment strategies for OS. However, further experiments are required to confirm these findings.

References