Identification of miRNA and genes involving in osteosarcoma by comprehensive analysis of microRNA and copy number variation data.

The aim of the present study was to understand the molecular mechanisms of osteosarcoma by comprehensive analysis of microRNA (miRNA/miR) and copy number variation (CNV) microarray data. Microarray data (GSE65071 and GSE33153) were downloaded from the Gene Expression Omnibus. In GSE65071, differentially expressed miRNAs between the osteosarcoma and control groups were calculated by the Limma package. Target genes of differentially expressed miRNAs were identified by the starBase database. For GSE33153, PennCNV software was used to perform the copy number variation (CNV) analysis. Overlapping of the genes in CNV regions and the target genes of differentially expressed miRNAs were used to construct miRNA-gene regulatory network using the starBase database. A total of 149 differentially expressed miRNAs, including 13 downregulated and 136 upregulated, were identified. In the GSE33153 dataset, 987 CNV regions involving in 3,635 genes were identified. In total, 761 overlapping genes in 987 CNV regions and in the genes in 7,313 miRNA-gene pairs were obtained. miRNAs (hsa-miR-27a-3p, hsa-miR-124-3p, hsa-miR-9-5p, hsa-miR-182-5p, hsa-miR-26a-5p) and the genes [Fibroblast growth factor receptor substrate 2 (FRS2), coronin 1C (CORO1C), forkhead box P1 (FOXP1), cytoplasmic polyadenylation element binding protein 4 (CPEB4) and glucocorticoid induced 1 (GLCCI1)] with the highest degrees of association with osteosarcoma development were identified. Hsa-miR-27a-3p, hsa-miR-9-5p, hsa-miR-182-5p, FRS2, CORO1C, FOXP1 and CPEB4 may be involved in osteosarcoma pathogenesis, and development.

Introduction

Osteosarcoma is the most common type of primary bone tumor in adolescents and young adults (1,2), and is characterized by an abundance of genomic aberrations (3). According to previous studies, the 5-year survival rate of osteosarcoma is 60–70%, and the prognosis has not significantly improved over the last 30 years (4). Despite improvements in osteosarcoma treatment, the molecular mechanism underlying osteosarcoma development remains unclear. Therefore, it is important to explore the molecular mechanism of osteosarcoma development to additionally improve osteosarcoma treatment. However, the mechanisms of osteosarcoma development are complex, and numerous factors, including genes associated with osteosarcoma (3), copy number variations (CNVs) across the whole genome (5) and miRNAs (6), may contribute to the development of osteosarcoma. In particular, CNV may contribute to the development of various types of cancer, including osteosarcoma and B-cell lymphoma (7,8). CNV has been considered as a marker for cancer predisposition (9). CNV may contribute to the pathogenesis of osteosarcoma (10) and the disease risk for diffuse large B-cell lymphoma (11).

Furthermore, a number of previous studies revealed that microRNAs (miRNAs/miRs) are associated with the development of osteosarcoma (12,13), and may be used as the molecular targets for osteosarcoma, including miR-143, miR-382 and miR-214. Downregulation of miR-143 is associated with the lung metastasis of osteosarcoma via the upregulated expression of matrix metalloproteinase 13 (14). Overexpression of miR-382 may suppress the metastasis of osteosarcoma (15). Upregulated expression of miR-214 may contribute to the pathogenesis of osteosarcoma and may be associated with adverse prognosis (16). In addition, miRNAs, including miR-143 (17), miR-27a (18), miR-223 (19), miR-191 (20), miR-133a (21) and miR-26a (22), may also be associated with the development of osteosarcoma. These results indicated that miRNAs may serve an important role in the pathogenesis of osteosarcoma.

Previously, miRNAs expression profiles (23,24) and SNP microarray (1) were used to identify the genes, miRNAs, and single nucleotide polymorphisms (SNPs) associated with the development of osteosarcoma. In certain studies, specific genes associated with osteosarcoma development were identified by the integrative analysis of copy number and gene expression data (25). Based on the aforementioned, comprehensive analysis of microRNA (GSE65071) and CNV (GSE33153) microarray data were performed in the present study to identify the miRNAs, and genes involved in the pathogenesis of osteosarcoma. In addition, a miRNA-gene regulatory network of the overlapping genes of the genes in CNV regions and target genes of differentially expressed miRNAs was constructed. These miRNAs and genes associated with osteosarcoma development were identified in order to additionally explore the mechanism of osteosarcoma, and may be used as candidate target miRNAs and genes for gene therapy.

Materials and methods

Microarray data

Microarray data (GSE65071 and GSE33153) were downloaded from Gene Expression Omnibus (GEO, http://www.ncbi.nlm.nih.gov/geo/). For GSE65071, there were 20 osteosarcoma samples and 10 healthy control samples, and the platform was GPL19631Exiqon human V3 microRNA PCR panel I+II. A total of 32 osteosarcoma samples were available on GSE33153, of which the platform was GPL6801 [GenomeWideSNP_6] Affymetrix Genome-Wide Human SNP 6.0 Array.

Data preprocessing and analysis of copy number variation

GSE65071 was preprocessed using preprocessCore package v1.38.1 (http://www.bioconductor.org/packages/release/bioc/html/preprocessCore.html) in R language. Then, the differentially expressed miRNAs between the osteosarcoma and control groups were calculated using the Limma package v3.32.2 (http://www.bioconductor.org/packages/release/bioc/html/limma.html), and the adjusted P-value <0.05 and |log (fold change) |>2 were chosen as the cut-off criterion. Then, two-way clustering analysis was performed using the gplots package v3.0.1 (https://cran.r-project.org/web/packages/gplots/index.html).

For GSE33153, the log R Ratio (LRR) and B allele frequency (BAF) were extracted using the affy2sv package v1.0.12 (https://bitbucket.org/brge/affy2sv), and then CNV calling was performed using PennCNV v2014 (07) (http://penncnv.openbioinformatics.org/). CNV identified in ≥4 samples was considered as CNV regions, which may be involved in the development of osteosarcoma.

Identification of target genes of miRNAs

The regulatory associations between long non-coding RNA, miRNA, competing endogenous RNA, mRNA and RNA binding proteins may be identified using the starBase database (http://starbase.sysu.edu.cn/), which included the clinical data and expression profiles of 14 types of cancer. The target genes of differentially expressed miRNAs were identified using the starBase database. The cut-off criterion of miRNA-gene pairs, which were used to screen the target genes, was listed as follows: i) miRNA-gene pairs were confirmed by ≥2 experiments; and ii) the change in the trend of miRNA and target gene expression values was the opposite.

Comprehensive analysis of miRNAs and CNV regions

Genes in the CNV regions were screened using scan_region.pl of PennCNV. The overlap between genes contained within the CNV regions and targets of the differentially expressed miRNAs were obtained. miRNA-gene pairs that overlapped were obtained based on the starBase database, and a miRNA-gene regulatory network was constructed.

Results

Analysis of differentially expressed miRNA and CNV regions

A total of 149 differentially expressed miRNAs between the osteosarcoma group and the control group were identified, including 13 downregulated, and 136 upregulated differentially expressed miRNAs (Table I). The clustering plot indicated that the differentially expressed miRNAs evidently separated the osteosarcoma samples from the control samples (Fig. 1).

Table I.

Top ten differentially expressed miRNAs.

miRNAs	logFC	P-value	Adjusted P-value
hsa-miR-624-5p	4.998135150	2.08×10−30	7.41×10−28
hsa-miR-505-5p	3.631886325	2.96×10−28	5.27×10−26
hsa-let-7f-2-3p	3.613251767	4.96×10−28	5.89×10−26
hsa-miR-877-3p	4.267538662	4.38×10−27	3.89×10−25
hsa-miR-183-5p	4.469024403	1.94×10−26	1.26×10−24
hsa-miR-342-5p	5.006618179	2.12×10−26	1.26×10−24
hsa-let-7f-1-3p	3.182939079	8.79×10−26	4.47×10−24
hsa-miR-671-5p	5.476974410	1.58×10−25	7.02×10−24
hsa-miR-95	5.163185811	2.96×10−25	1.13×10−23
hsa-miR-499a-5p	5.429803038	3.18×10−25	1.13×10−23

miRNAs/miR, microRNA; hsa, Homo sapiens.

Figure 1.

Clustering plot of the differentially expressed microRNAs when comparing the osteosarcoma group and the control group. Green and red colors represent the low, and high expression values, respectively.

For GSE33153, 987 CNV regions were identified. The distribution of the 987 CNV regions among chromosomes is presented in Fig. 2. In Fig. 2, the red and blue column represents the deletions, and overlap of copy number, respectively. The longer the column, the higher the number of samples that exhibited deletions or overlapping in these CNV regions.

Figure 2.

Distribution of 987 CNV regions among chromosomes. Red and blue column represented the deletions, and overlap of copy numbers, respectively. The longer column, the higher the number of samples that exhibited deletions or overlapping in these CNV regions. CNV, copy number variation.

Analysis of target genes of differentially expressed miRNAs

In total, 7,313 miRNA-gene pairs were screened for 149 differentially expressed miRNAs, and 5 differentially expressed miRNAs (hsa-miR-200b-3p, hsa-miR-200a-3p, hsa-miR-429, hsa-miR-34a-5p, and hsa-miR-9-5p) with corresponding target genes were identified, as summarized in Table II.

Table II.

miRNA-target gene pairs.

miRNAs	Counts	Target genes[a]
hsa-miR-200b-3p	282	ICK, NPTX1, DPY19L1, RECK, TRIO, DDX3Y, GALNT2, R3HDM2, PPM1B, TIMP2
hsa-miR-200a-3p	159	GIGYF1, GIGYF1, MTF2, CASC4, ULK2, MBNL3, KPNA4, ZEB2, ZEB2, TET3
hsa-miR-429	235	SEC24A, TRIM52, DPY19L1, PAPOLA, GALNT2, KIAA1432, MEX3C, PLCG1, LRIG1, DCBLD2
hsa-miR-34a-5p	92	TNRC18, SIDT2, ARPP19, NAA50, RCAN1, DNAJB1, XPO5, HCN3, SIRT1, MBD6
hsa-miR-9-5p	144	ULK2, SMARCE1, ZC3H10, VAT1, BCAT2, DICER1, ZC3H12A, GTPBP3, SH2B3, ITM2B

The genes listed do not cover the exhaustive group of target genes for each miRNA. miRNAs/miR, microRNA; hsa, Homo sapiens. Counts represent the number of target genes which were contained within the CNV regions.

Comprehensive analysis of differentially expressed miRNA and CNV regions

A total of 761/3,635 genes in 987 CNV regions and the genes of 7,313 miRNA-gene pairs were identified overlap. The miRNAs with highest degrees of overlap (hsa-miR-27a-3p degree, 112; hsa-miR-124-3p degree, 102; hsa-miR-9-5p degree, 90; hsa-miR-182-5p degree, 79; hsa-miR-26a-5p degree, 79) and genes (FRS2 degree, 14; CORO1C degree, 12; FOXP1 degree, 11; CPEB4 degree, 11; GLCCI1 degree, 10) are presented in Table III.

Table III.

miRNAs and genes in the miRNA-gene regulatory network.

A, miRNAs in the miRNA-gene regulatory network
Marker	Degrees	Counts	Target genes[a]
hsa-miR-27a-3p	112	51	TNRC18, PPARA, ACLY, DNAJB9, NR2F2, HMGCS1, SGMS1, CKAP4, RPN2, CALM
hsa-miR-124-3p	102	42	SGMS1, DRAM1, SUCLG2, FRS2, NUDCD2, FAR1, PIP4K2C, SERTAD3, SH3PXD2A, KLHL24
hsa-miR-9-5p	90	39	FREM2, SMARCE1, KLHDC10, DYRK1B, SLC39A14, FURIN, BAHD1, MAPKAPK2, AP2M1, BCL6
hsa-miR-182-5p	79	38	RTN4, STARD13, FRS2, KIAA1217, CNOT6, EXOC4, SNAP23, KCMF1, QSER1, SYPL1, BDNF
hsa-miR-26a-5p	79	29	SRP19, MARK1, SEMA6D, ACSL3, EIF4G2, LSM12, MAPK6, CCNJL, MFHAS1, COX5A
hsa-miR-429	79	32	SUZ12, CRKL, TOB1, FXR1, LIN7B, EVI5, GLCCI1, FRS2, ARID4B, SSH2
hsa-miR-141-3p	71	30	RHEB, LENG8, ATP2A2, BAHD1, GLCCI1, ATP1B1, PTPRG, MBTD1, GRIN2D, STAT5B
hsa-miR-96-5p	70	33	ARPP19, CELSR1, KIAA1217, DOCK1, CAPNS1, CCNG1, PROK2, APPL1, PGAP1, SHC1
hsa-miR-34a-5p	63	30	HCN3, HECW2, STC1, MAP2K1, APH1A, NUMBL, NFE2L1, GREM2, ARID4B, LGR4
hsa-miR-200a-3p	63	27	RHEB, NCKAP5, B3GNT5, LENG8, CALU, ATP2A2, BAHD1, NRP1, SPAG9, GLCCI1
B, Genes in the miRNA-gene regulatory network
Marker	Degrees	Counts	Target genes
FRS2	14	–	–
CORO1C	12	–	–
FOXP1	11	–	–
CPEB4	11	–	–
GLCCI1	10	–	–
CELF1	10	–	–
MET	10	–	–
ZFHX4	9	–	–
DOCK4	8	–	–
MYH10	8	–	–

The genes listed do not cover the exhaustive group of target genes for each miRNA. miRNAs/miR, microRNA; hsa, Homo sapiens. Counts represent the number of target genes which were contained within the CNV regions.

Discussion

Osteosarcoma is characterized by an abundance of genomic aberrations (3). In the present study, a comprehensive analysis of microRNA data and CNV microarray data was performed to identify miRNAs, and genes associated with the pathogenesis of osteosarcoma. A total of 149 differentially expressed miRNAs, including 13 down- and 136 upregulated, in the GEO GSE65071 dataset were identified. For the GEO GSE33153 dataset, 987 CNV regions were identified. In addition, a miRNA-gene regulatory network of 761 overlapping genes out of 3,635 genes in 987 CNV regions and the genes in 7,313 miRNA-gene pairs was constructed. Concurrently, miRNAs (hsa-miR-27a-3p, hsa-miR-124-3p, hsa-miR-9-5p, hsa-miR-182-5p, hsa-miR-26a-5p) and genes [fibroblast growth factor receptor substrate 2 (FRS2), coronin 1C (CORO1C), forkhead Box P1 (FOXP1), cytoplasmic polyadenylation element binding protein 4 (CPEB4) and glucocorticoid induced 1 (GLCCI1)] with the highest degree of overlap in the miRNA-gene regulatory network were identified. Higher degrees represent an increased level of association to the development of osteosarcoma. Therefore, hsa-miR-27a-3p, hsa-miR-124-3p, hsa-miR-9-5p, hsa-miR-182-5p, hsa-miR-26a-5p, FRS2, CORO1C, FOXP1, CPEB4 and GLCCI1 may contribute to the pathogenesis of osteosarcoma.

The expression of hsa-miR-27a-3p serves an important role in the development of a number of types of cancer, including breast and pancreatic cancer. The expression of hsa-miR-27a-3p was positively correlated with a hypoxia gene signature in breast cancer (26). The inhibition of hsa-miR-27a-3p expression may exhibit potentially antiproliferative effects in pancreatic cancer (27). In addition, hsa-miR-27a-3p was considered as a candidate biomarker for Alzheimer's disease (28). However, the correlation between hsa-miR-27a-3p and osteosarcoma has not been widely investigated. In the present study, hsa-miR-27a-3p was identified as a hub gene with a high degree of association in miRNA-gene regulatory network. In combination with previous studies, the present study hypothesized that hsa-miR-27a-3p may serve a role in the development of osteosarcoma.

Hsa-miR-9-5p, hsa-miR-155-5p and hsa-miR-203 are potent prognostic factors for acute myeloid leukemia (29). In addition, a 3-miRNA scoring system (hsa-miR-9-5p, hsa-miR-155-5p and hsa-miR-203) was used for the prognostication of the patients with de novo acute myeloid leukemia (30). Hsa-miR-182-5p has been considered as a marker for distinguishing between human ovarian cancer tissues and normal tissues (31). The knockdown of miR-182-5p significantly decreased the growth of prostate tumor, and FOXF2, reversion inducing cysteine rich protein with kazal motifs (RECK) and metastasis suppressor 1 were identified as potential target genes of miR-182-5p (32). Previous studies have also revealed that mothers against decapentaplegic homolog 4 (Smad4) and RECK were the potential target genes of miR-182-5p, and that miR-182-5p served a key role in bladder cancer by knocking down the expression of Smad4 and RECK (33). In the present study, hsa-miR-9-5p and hsa-miR-182-5p were identified as the hub genes with high degrees of association in the miRNA-gene regulatory network, implying that hsa-miR-9-5p and hsa-miR-182-5p may serve a key role in the pathogenesis of osteosarcoma.

In the miRNA-gene regulatory network, FRS2, CORO1C, FOXP1, CPEB4 and GLCCI1 with high degrees of association were identified, implying that these genes may be associated with the regulatory mechanism of osteosarcoma. FRS2 is considered to be a gene that is associated with numerous types of cancer, including ovarian cancer, liposarcoma and prostate cancer (34–36). FRS2 serves as an amplified oncogene that induces the downstream activation of the Ras-mitogen-activated protein kinase pathway in high-grade serous ovarian cancer (36). In addition, FRS2 serves an essential role in fibroblast growth factor receptor (FGFR) signaling, and activated FGFR/FRS2 signaling may lead to the development of high-grade liposarcoma (34). Concomitantly, a previous study also identified that the functional overlap of FRS2 and FRS3 may mediate mitogenic FGF signaling in prostate cancer (37). From these data, the present study hypothesized that FRS2 was involved in the regulatory mechanism of osteosarcoma.

CORO1C is a target gene of miR-206 that has been demonstrated to inhibit cell migration in triple-negative breast cancer (38). CORO1C is also considered as a target gene of miR-1/133a following genome-wide gene expression and luciferase reporter assay analyses (39). Furthermore miR-1 and miR-133a inhibited the proliferation, migration, and invasion of lung-squamous cell carcinoma cells (39). In addition, miR-1 and miR-133b mediated cell proliferation, and cell cycle progression by regulating the expression of hepatocyte growth factor receptor protein in human osteosarcoma (40). miR-133b and miR-206 expression were identified to be significantly decreased and may be used as potential prognostic markers for patients with osteosarcoma (41). Decreased miR-206 expression is associated with the development of osteosarcoma, and the transfection of miR-206 mimics promoted cell apoptosis, and inhibited cell invasion and migration (42). CORO1C may serve as the target gene of miR-1, miR-133b and miR-206, which have all been associated with osteosarcoma development. Furthermore, CORO1C was identified as a node gene with a high degree of association in the miRNA-gene regulatory network. However, CORO1C was not identified as a target gene of the miRNAs in the present study. Together, the results indicate that CORO1C may serve a role in osteosarcoma pathogenesis and development by regulating miRNAs.

FOXP1 expression is associated with the development of osteosarcoma (43,44). CPEB4 serves a role in metastatic cancer and cancer progression (45,46); however, whether CPEB4 is involved in the pathogenesis of osteosarcoma remains unclear. In the present study, FOXP1 and CPEB4 were also identified as the node genes with high degrees of association in the miRNA-gene regulatory network. Based on the results, it was hypothesized that FOXP1 and CPEB4 may contribute to osteosarcoma progression, and development. However, these results require additional confirmation.

In the present study, hsa-miR-27a-3p, hsa-miR-9-5p, hsa-miR-182-5p, FRS2, CORO1C, FOXP1 and CPEB4 were identified as node genes with high degrees of association in the miRNA-gene regulatory network, and may serve a role in the pathogenesis and development of osteosarcoma.