Integrated bioinformatics analysis of miRNA expression in osteosarcoma

Abstract

The expression of miRNA influencing the pathogenesis of OS have been reported previously, however, different samples selection and sequencing platforms made obvious differences in miRNA expression analysis. We aim to identify reliable prognostic and treatment biomarkers for OS by systematic analysis of miRNAs expression data sets from biased data set. Seven miRNA data sets were selected from corresponding articles. Collectively, two miRNAs, hsa-miR-19-3p and hsa-miR106b-3p, and transcription factor SIX3 were identified and may be reliable markers for prognostic and treatment of osteosarcoma.

Keywords:

• miRNA expression
• osteosarcoma
• SIX3
• target genes

Introduction

Osteosarcoma (OS) remained the most common and destructive primary bone tumor that affected children and adolescent. It was characterized by malignant osteoid production and osteoblastic differentiation. Although OS have been treated with chemotherapy for decades, patients with recurrent or metastatic OS still had very poor prognosis (Ottaviani and Jaffe 2009, Tan et al. 2009).

To date, however, the molecular pathogenesis and etiology of OS were still not clearly elucidated. Micro-RNA (miRNA), a short conserved non-coding RNA with 22–24 nucleotides long, was believed to be a promising diagnostic and prognostic tool of malignant neoplasm (Doench and Sharp 2004). MiRNA bound to the sufficient target mRNA, which led to mRNA degradation or preventing mRNA from being translated. Therefore, miRNAs regulated target genes expression at post-transcriptional level. Over-expression of miRNAs usually caused the decreased expression of target genes. Evidence showed that 98 miRNAs located at fragile sites and genomic regions involved in cancer, indicating that miRNAs were closely related to tumorigenesis and development (Calin et al. 2004). MiRNAs played a critical role in tumor growth, progression, metastasis, and drug resistance through regulating the expression of oncogenes or tumor suppressor genes (Shi et al. 2008, Skaftnesmo et al. 2007). However, different sample selection and sequencing platform made an obvious discrepancy in miRNA expression analysis.

In this study, meta-analysis was used to determine the key miRNA in the process of the onset of OS. Then, target genes corresponding to the differential expression miRNAs (DEMs) were predicted and its functions were annotated. Furthermore, transcription factor analysis was analyzed using Tfacts database. This study demonstrated vital significance of miRNA in prevention and treatment of OS.

Materials and methods

Screening of miRNA expression data sets and DEMs in osteosarcoma

“Oseteosarcoma” and “miRNA” were used as key words to retrieve the existing articles of miRNA expression in OS in the Google academic web site (https://scholar.google.com), and the expression dates were carefully screened. Only articles in which miRNAs sequence of samples with OS and contrast were in detail were selected and subsequently analyzed. In these articles, specific DEMs were extracted, and statistics of up- and down-regulated DEMs compared to the control samples were listed. Number of DEMs was counted in independent data set. Only DEMs appearing at least three separate data sets were reserved. The overall workflow for this study is showed in Figure 1.

Figure 1. The analysis process of this study.

Target genes of DEMs in osteosarcoma

The sequence information of DEMs were retrieved in miRBase (Griffiths-Jones et al. 2008, Kozomara and Griffiths-Jones 2011) and collated in fasta files. Target Scan (Lewis et al. 2005) software was adopted to predict corresponding target genes of DEMs using the default parameters.

Functional enrichment analysis of the target genes of DEMs in osteosarcoma

Gene Ontology (GO) and Pathway analyses were used to annotate the forecasted DEMs by using DAVID database (Dennis et al. 2003). Top ten GO terms of biological processes in GO annotation results were extracted out, and pathways were screened with FDR value less than 0.05.

Transcription factor analysis of up- and down-regulated DEMs in osteosarcoma

Target genes corresponding to up- and down-regulated miRNA were analyzed by using Tfacts database (Essaghir et al. 2010) (http://www.tfacts.org/). Only transcription factors predicted with four indexes (P values, q-value, E-value and FDR) less than 0.05 were considered to be reliable transcription factors. Statistics on transcription factors of target genes of up- and down-regulated DEMs were made to select the shared and specific transcription factors between up- and down-regulated DEMs.

Results

Screening of differentially expressed miRNA

After retrieval, seven independent miRNA expression data sets were selected from the existing researched of miRNA expression in OS, in which the data set provided detailed miRNA expression date and up- and down-regulated DEMs compared to the control sample. For convenience, these seven miRNA expression data sets were numbered in detail as follows: data set 1 (Maire et al. 2011), data set 2 (Namløs et al. 2012), data set 3 (Jones et al. 2012), data set 4 (Duan et al. 2011), data set 5 (Hu et al. 2012), data set 6 (Zhao et al. 2013), and data set 7 (Thayanithy et al. 2012).

Detailed information including the type, size, the number of sample and the type and number of chip were listed (Table 1). Furthermore, DEMs were validated by qPCR in all articles of these seven data sets, indicating the results more reliable.

Table 1. Information of all seven data sets.

Data set	Samples	Assay type	Number of miRNA probes	Validated
1	Seven osteosarcoma tumor samples(OS4, OS7, OS9, OS11, OS15, OS17, OS18)	Agilent	799	qRT-PCR
2	19 osteosarcoma cell lines	Agilent	799	qRT-PR
3	OS versus control tissues	Agilent	1100	qRT-PR
4	osteosarcoma cell lines	LC Sciences	–	qRT-PR
5	Osteosarcoma cell line (MG – 63) and osteoblast cell line	Illumina	1146	qRT-PCR
6	23 patients with osteosarcoma, 3 normal muscle tissues	miRCURY LNA	–	qRT-PCR
7	OS versus control tissues	Illumina	1135	qRT-PCR

Distribution of DEMs in all seven data sets were analyzed using SVG. The results showed that there were big differences in both up- and down-regulated miRNAs between different data sets (Figure 2). Statistics were further made to demonstrate the number of up- and down-regulated miRNA in each data set (Figure 3).

Figure 2. Distribution of DEMs alterations in OS as reported by primary studies. Short vertical bars indicate down- or up-regulated miRNA. There were significant differences in both up-regulated and down-regulated miRNAs between different data sets.

Figure 3. Statistics of the number of up- and down-regulated miRNA in the existing seven independent miRNA expression data set in osteosarcoma the left pillars in each dataset represents the quantity of up-regulated miRNA compared with the normal sample tissue expression; the right pillars in each dataset represents the quantity of down-regulated miRNA.

In the seven existing OS miRNA expression data sets, numbers of DEMs were different. Numbers of DEMs in five data sets were more than 30, however, numbers of DEMs in other two data sets were less than 10, that were data set 4 with 5 DEMs and data set 5 with 9 DEMs, respectively. In these seven existing miRNA expression data sets, up-regulated DEMs were most in data set 2 with 40 up-regulated DEMs, followed by the data set 7 with 22 up-regulated DEMs. Data set 2 possessed the least up-regulated DEMs, only 2. The other hand, data set 1 and data set 6 possessed the same and largest number of down-regulated DEMs, 29, followed by data set 3 with 27. In data set 4, the number of down-regulated DEMs was the least, only 3.

Screening and location of DEMs

DEMs were extracted from the corresponding literature compared with the control samples (Figure 4). In these seven independent data sets, there were 65 non-redundant up-regulated miRNAs and 76 non-redundant down-regulated miRNAs, which corresponded to different data sets. To up-regulated miRNAs, total number of data set specific up-regulated miRNAs was 40, accounting for 61.5% of the total number of non-redundant up-regulated miRNAs. On the other hand, total number of data set specific down-regulated miRNAs was 51, accounting for 67.1% of total down-regulated non-redundant miRNAs. No matter up-regulated or down-regulated miRNA, proportion of data set specific miRNA in non-redundant miRNAs were more than 60%, which indicated that the existing sequencing results in OS miRNA expression profile showed great difference from each other. Therefore, it was significant to analyze the miRNA expression in OS by bioinformatics analysis.

Figure 4. Number of up- and down-regulated miRNA correspond to the seven independent data set.

Based on the support of data sets, reliable DEMs were selected, which were supported by 3 data sets at least. According to this standard, 4 up-regulated and 8 down-regulated miRNAs were screened to be high trustworthy miRNAs. Detailed location information of miRNAs was analyzed by using miRBase database (Table 2). Of these 12 DEMs, there was 3 miRNAs in chromosome 14 with the largest number of miRNAs, accounting for 25% of DEMs, followed by chromosome 7, chromosome 13 and chromosome 19, with two miRNA. There was one miRNA in chromosome 9, chromosome 17 and chromosome 18, respectively.

Table 2. List of miRNA with meta characteristics in osteosarcoma.

miRNA	chromosome	beg	end	strand	supporteddata sets
Upregulated
hsa-miR-106	chr7	100093993	100094074	−	3
hsa-miR-17	chr13	91350605	91350688	+	3
hsa-miR-181	chr19	13874699	13874808	+	8
hsa-miR-19	chr13	91350891	91350972	+	3
Downregulated
hsa-miR-126	chr9	136670602	136670686	+	4
hsa-miR-127	chr14	100882979	100883075	+	3
hsa-miR-133	chr18	21825698	21825785	−	4
hsa-miR-150	chr19	49500785	49500868	−	3
hsa-miR-335	chr7	130496111	130496204	+	4
hsa-miR-376	chr14	101040782	101040849	+	7
hsa-miR-409	chr14	101065300	101065378	+	3
hsa-miR-451	chr17	28861369	28861440	−	3

Target gene prediction of DEMs

Target genes of these 12 DEMs were predicted using TargetScan software (Cambridge, MA). Though all parameters used were default, the numbers of predicted target genes were different. Number of predicted target genes by the hsa-miR-19–3p was the most with 1330, followed by hsa-miR106b-3p with 698 target genes. According to the default parameters, the numbers of predicted target genes of the other 10 DEMs were 100. According to the number of predicted target genes, miRNA possessed different conservative. As we all know, the more miRNA were conservative, the more target genes were in the genome. According to the results, hsa-miR-19–3p and hsa-miR106b-3p were highly conservative, but the other 10 miRNAs had the same level of conservatism.

Function analysis of target genes of DEMs

Enriched metabolic pathways of target genes regulated by each specific miRNAs were extracted out (Figure 5). Pathways enrichment of target genes corresponding to seven of twelve DEMs were established and enrichment level of metabolic pathways were different. According to the analysis, target genes of miR-106-b-3p possessed the most extensive pathway enrichment and the highest pathway enrichment degree.

Figure 5. Metabolic pathways enrichment of target genes of miRNA Target genes of Seven out of 12 DEMs possessed appropriate Pathway enrichment. The gradient from green to red represents the corresponding Pathway proportion (%) of the terms in the miRNA target genes, and the deeper the color, the higher the proportion.

Functional annotations of target genes of up- and down-regulated miRNA were analyzed by DAVID. Top ten of significant difference biological processes and pathways with FDR <0.05 were selected. Significantly enriched processes of up-regulated miRNA were associated with regulating transcription and RNA metabolism in OS (Table 3), indicating that target genes of up-regulated miRNA were primarily transcription factor by strengthening the regulation of transcription factor to enhance regulation of gene expression. Target genes of down-regulated miRNA played a role in regulating phosphorus metabolic processes and secretion of tumor necrosis factor, however, the biological significance of them in OS need further study (Table 4).

Table 3. GO annotation of target genes of up-regulated miRNA in osteosarcoma.

Category	Term	Count	%	P value	Foldenrichment	FDR	Fisher exact
GOTERM_BP_FAT	regulation of transcription	438	2.2	4.90E-25	1.5	9.00E-22	3.00E-25
GOTERM_BP_FAT	transcription	360	1.8	2.40E-21	1.6	4.50E-18	1.40E-21
GOTERM_BP_FAT	regulation of RNA metabolic process	318	1.6	3.10E-20	1.6	5.70E-17	1.80E-20
GOTERM_BP_FAT	regulation of transcription, DNA-dependent	312	1.6	4.70E-20	1.6	8.70E-17	2.70E-20
GOTERM_BP_FAT	regulation of transcription from RNA polymerase II promoter	141	0.7	5.50E-12	1.8	1.00E-08	2.70E-12
GOTERM_BP_FAT	homophilic cell adhesion	42	0.2	1.90E-10	2.9	3.40E-07	4.90E-11
GOTERM_BP_FAT	positive regulation of transcription, DNA-dependent	95	0.5	6.10E-09	1.8	1.10E-05	3.00E-09
GOTERM_BP_FAT	positive regulation of RNA metabolic process	95	0.5	9.50E-09	1.8	1.80E-05	4.70E-09
GOTERM_BP_FAT	positive regulation of gene expression	108	0.5	2.50E-08	1.7	4.60E-05	1.30E-08
GOTERM_BP_FAT	positive regulation of nucleobase, nucleoside, nucleotide and nucleic acid metabolic process	114	0.6	2.60E-08	1.7	4.80E-05	1.40E-08

Table 4. GO annotation of target genes of down-regulated miRNA in osteosarcoma.

Category	Term	Count	%	P value	Foldenrichment	FDR	Fisher exact
GOTERM_BP_FAT	regulation of phosphorus metabolic process	36	0.6	1.20E-04	2	2.10E-01	5.80E-05
GOTERM_BP_FAT	regulation of phosphate metabolic process	36	0.6	1.20E-04	2	2.10E-01	5.80E-05
GOTERM_BP_FAT	response to tumor necrosis factor	5	0.1	1.00E-03	10.4	1.80E + 00	6.90E-05
GOTERM_BP_FAT	positive regulation of cell proliferation	32	0.5	1.60E-04	2.1	2.70E-01	7.00E-05
GOTERM_BP_FAT	regulation of phosphorylation	34	0.5	2.70E-04	2	4.70E-01	1.30E-04
GOTERM_BP_FAT	regulation of kinase activity	28	0.5	3.50E-04	2.1	6.10E-01	1.50E-04
GOTERM_BP_FAT	regulation of protein kinase activity	27	0.4	4.70E-04	2.1	8.30E-01	2.10E-04
GOTERM_BP_FAT	regulation of transferase activity	28	0.5	6.60E-04	2	1.20E + 00	3.10E-04
GOTERM_BP_FAT	positive regulation of fatty acid beta-oxidation	3	0	1.30E-02	16.2	2.00E + 01	4.80E-04
GOTERM_BP_FAT	protein heterooligomerization	8	0.1	2.90E-03	4.2	4.90E + 00	6.00E-04

Enriched metabolic pathways of target genes of miRNA were selected with FDR <0.05 (Figure 6). The results showed that target genes of up-regulated miRNA were mainly associated with six main metabolic pathways, namely endocytosis, MAPK signaling pathways, glioma, long-term potentiation, ErbB signaling pathways and cancer signaling pathways, and the number of genes in cancer signaling pathway enrichment was the largest, suggesting that up-regulated miRNA influenced the pathogenesis of OS through participating in a variety of metabolic pathways to regulate the expression of target genes. Analysis showed that there was no significantly enriched pathway of target genes corresponding to down-regulated miRNA in OS.

Figure 6. Metabolic pathways enrichment analysis of target genes of up-regulated miRNA in osteosarcoma Target genes of up-regulated miRNA were mainly enriched in four pathways, with cancer possessing the highest enrichment degree.

Transcription factors analysis of target genes of DEMs

Transcription factors corresponding to target genes of up- and down-regulated miRNA were analyzed, and contrastive analysis was carried out between them. As the results showed, there were 136 transcription factors corresponding to target genes and 622 interactions between transcription factors and target genes in up-regulated miRNA analysis. On the other hand, 81 transcription factors corresponding to target genes and 622 interactions between transcription factors and target genes were involved in down-regulated miRNA (Figure 7(A)). Contrastive analysis showed that of the total 152 transcription factors, 65 were shared in the two types of target genes, accounted for 42.8% of the total transcription factors. Furthermore, transcription factors corresponding to target genes of up-regulated miRNA showed higher specificity, consistent with the results of GO annotations, in which target genes of up-regulated miRNA enriched in regulating transcription. According to the analysis, main transcription factors with credibility E-value and cross ratio were not greater than 0.05, only transcription factor SIX3 had high cross ratio, indicating that the identification of transcription factors were reliable (Figure 7(B)).

Figure 7. Transcription factor analysis of target genes of up- and down-regulated miRNAs. (A) Share and specific transcription factor comparison of target genes of up- and down-regulated miRNA. The left background represents number of transcription factors of target genes of up-regulated miRNA. The right background on behalf of number of transcription factors of target genes of down-regulated miRNA. (B) E-value and cross ratio statistics of important regulating transcription factors.

Discussion

In this study, seven independent miRNA expression data sets were selected from the existing researches related to miRNA expression in OS and carried out with bioinformatics analysis. Then up- and down-regulated miRNAs were screened and analyzed. Two miRNA, hsa-miR-19–3p and hsa-miR106b-3p, were identified and might be a predicted signature of osteosarcoma. Target genes of DEMs were certificated, and the functions were analyzed. Finally, transcription factors corresponding to target genes were screened. This study showed a vital significance of miRNA in prevention and treatment of OS.

According to our results, 4 up-regulated and 8 down-regulated miRNAs were screened to be high trustworthy miRNAs. Two miRNAs, hsa-miR-19–3p and hsa-miR106b-3p, were selected as biological markers which might be potential signatures for prevention and treatment of OS. Hsa-miR-19a-3p was the oncogenic component of the mir-17/92 cluster (Olive et al. 2009, Ventura et al. 2008). When differential expression of microRNAs was detected in plasma of patients with laryngeal squamous cell carcinoma, hsa-miR-19a-3p was confirmed as a potential early-detection marker for laryngeal squamous cell carcinoma (Ayaz et al. 2013). Also, hsa-miR-19a-3p was texted in breast cancer grouped by clinicopathological subtypes (Liu et al. 2013). Down-regulation of miR-181b promoted apoptosis by targeting CYLD in thyroid papillary cancer, and hsa-miR-19a-3p was detected and played an important role (Li et al. 2014). Therefore, hsa-miR-19–3p might play a vital role in OS, and be a potential treatment marker. Hsa-miR106b-3p was only detected as the circulating miRNA in myasthenia gravis. In this study, a total of 698 target genes were highly associated with hsa-miR106b-3p, which indicated that hsa-miR106b-3p might be an important biomarker. However, the exact function in OS need further study.

Significantly enriched processes of up-regulated miRNAs were involved with regulating transcription and RNA metabolism in OS, indicating that target genes of up-regulated miRNAs were primarily transcription factor by strengthening the regulation of transcription factor to enhance regulation of gene expression. As we all know, up-regulation of gene transcription, especially genes related with cell proliferation and carcinogenesis, might lead to cell transform to cancer cell. Regulation of OS oncogenesis process by miRNA was a comprehensive feedback adjustment process. Target genes of down-regulated miRNA played a role in regulating main phosphorus metabolic processes and secretion of tumor necrosis factor the process. Compared to normal cells, choline kinase and phospholipase C were significantly overexpressed, and lysophospholipase 1, phospholipase A2, and phospholipase D were significantly under-expressed, in breast cancer cells (Glunde et al. 2004). The overall activity of phosphatidylcholine-specific phospholipase C and phospholipase D was also higher (∼5-fold) in human epithelial ovarian carcinoma cells (EOC), suggesting that both biosynthetic and catabolic pathways of the phosphatidylcholine cycle likely contributed to phosphocholine accumulation. Evidence of abnormal phosphatidylcholine metabolism might have implications in EOC biology and might provide an avenue to the development of noninvasive clinical tools for EOC diagnosis and treatment follow-up (Iorio et al. 2005). Collectively, abnormal phosphorus metabolic processes might contribute to the oncogenesis of OS. Researches had shown that cannabinoids induced cancer cell proliferation via tumor necrosis factor alpha-converting enzyme (TACE/ADAM17)-mediated transactivation of the epidermal growth factor receptor (Hart et al. 2004). Epidermal growth factor receptor transactivation mediated tumor necrosis factor-induced hepatocyte replication (Argast et al. 2004). Therefore, changes of tumor necrosis factor were highly associated with cell proliferation and thus influenced oncogenesis of OS. However, the biological significance of them in OS needed further analysis. Collectively, analyses of target genes further proved that the selected miRNAs were highly associated with the oncogenesis of OS and might be potential biomarker in treatment of OS.

Pathway enrichment analysis showed that target genes of up-regulated miRNA were mainly associated with six main metabolic pathways, namely endocytosis, MAPK signaling pathways, glioma, long-term potentiation, ErbB signaling pathways and cancer signaling pathways, which were consistent with the results of analysis of target genes. However, there was no significantly enriched pathway of target genes of down-regulated miRNA OS.

SIX3 was a human homologue of the highly conserved sine oculis gene family essential during embryonic development in vertebrates, and encoded a homeo-domain containing transcription factor (Mo et al. 2013). Researches were reported that the homeotic protein Six3 was a co-activator of the nuclear receptor NOR-1 and a co-repressor of the fusion protein EWS/NOR-1 in human extraskeletal myxoid chondrosarcomas (Laflamme et al. 2003). Wnt1 downstream signaling components had been well studied and activated in human cancer. Reports showed that MTA1s and MTA1 as important upstream modifiers of the Wnt1 transcription, and consequently its functions, by directly inhibiting the transcription of Six3, allowed derepression of Wnt1 transcription (Kumar et al. 2010). Study showed that SIX3 might play an important role as a novel suppressor in human lung cancer. SIX3 had potential as a novel prognostic biomarker for patients with lung adenocarcinomas (Mo et al. 2013). In our study, transcription factor SIX3 had high cross ratio, indicating that the identification of transcription factors were reliable. Collectively, SIX3 might be an important transcription factor in oncogenesis of OS and might be a vital biomarker for prevention and treatment of OS.

Collectively, two miRNAs, hsa-miR-19–3p and hsa-miR106b-3p were screened and might be reliable signatures for prevention and treatment of osteosarcoma. Target genes function annotation analyses showed that the up-regulated miRNA target genes were mainly involved with transcription factors, and down-regulated miRNA were related to phosphorus metabolism regulation and secretion of tumor necrosis factor. Transcription factor SIX3 with high cross ratio was screened and might also be reliable signatures of osteosarcoma. However, all this conclusions were deduced by bioinformatics analysis and need further study to confirm them.

Disclosure statement

All authors declared no conflict of interest.