Bioinformatics Analysis of Neuroblastoma miRNA Based on GEO Data.

OBJECTIVE: To analyze the changes in downstream genes, signaling pathways, and proteins based on the difference of microRNA (miRNA) expression in neuroblastoma (NB).
METHODS: GSE128004 second-generation sequencing expression data were downloaded from GEO, and Limma package of R language was used to analyze differential expression, and a volcano map and heat map were drawn; the target genes corresponding to the differential miRNA were found using the miWalk web tool, and GO (Gene Ontology) and KEGG (Kyoto Encyclopedia of Genes and Genomes) were performed. The key genes were identified and verified in the TCGA database.
RESULTS: A total of 34 differentially expressed miRNAs were screened out. Among them, 22 up-regulated miRNAs predicted 1163 target genes and 12 down-regulated miRNAs predicted 1474 target genes. Target genes were enriched and analyzed by KEGG to find the FOXO signal pathway, mTOR signal pathway, AMPK signal pathway, and other signal pathways. After GO analysis, axon formation, regulation of chemical synaptic transmitters, regulation of nerve synapses, regulation of cross-synaptic signals, and other physiological processes were assessed. A total of 16 key genes were obtained by PPI analysis, and the survival analysis of TP53 and ATM genes verified in the TCGA database showed statistical significance.
CONCLUSION: The 34 differential miRNAs may be related to the occurrence and development of NB. TP53 and ATM are related to the prognosis of NB. The role and mechanism of TP53 and ATM in NB need to be further verified.

Introduction

Neuroblastoma (NB) is the most common extracranial solid tumor in children. The incidence is about 25–50 cases per 100 million people. The median age of diagnosis is 18 months,1 and about 90% of tumors occur before the patient is 10 years old.2 The incidence in China is on the rise, with an average annual incidence of about 3000 cases.3 NB is a neuroendocrine tumor that occurs in the developing sympathetic nervous system, from any neural crest element, resulting in adrenal or sympathetic ganglion tumors or both. Although a few NBs may subside without treatment, most cases have no specific clinical manifestations, and metastasis has often already occurred by the time of diagnosis, leaving the patient with poor prognosis. Early diagnosis and treatment are very important. Although several advanced investigations that include next-generation sequencing of tumor cells and immunological assessments to define targets in the tumor microenvironment can contribute to the diagnosis and staging of NB, these detective methods are invasive, time-consuming and also difficult to repeat the procedure.4 More recently, there has been a growing number of publications focusing on the function of exosomes in body fluids as an early and non-invasive diagnostic biomarker in malignant tumors. For example, serum exosomal hsa-miR-18a, hsa-miR-221, hsa-miR-222, and hsa-miR-224 have been reported as novel biomarkers for hepatocellular carcinoma (HCC),5 and these serum exosomal miRNAs are better in distinguishing HCC from chronic hepatitis B or liver cirrhosis compared to the corresponding serum circulating miRNAs. However, there are currently few studies on miRNAs in neuroblastoma. In this study, we analyzed the second-generation sequencing data of patients with NB in the public database GEO to identify differentially expressed microRNAs (miRNAs) and their target genes in the pathogenesis of NB.

Materials and Methods

Data Collection

The original data of microRNA second-generation sequencing were collected from the GEO database (). The ID number of the second-generation sequencing data is GSE128004, and the sequencing platform was GPL23227 BGISEQ-500 (Homo sapiens). The sequencing samples were collected from the peripheral blood of 15 children with NB and 3 healthy controls, including 15 NB samples (GSM3660254-GSM3660268) and 3 control samples (GSM3660269-GSM3660271).

Data Preprocessing and Differential Expression Analysis

After downloading the processed sequencing expression matrix, the Limma package6 was used to analyze the differences in miRNA expression between normal samples and NB samples, and a volcano map and heat map were drawn. Genes with incomplete annotation information, very low expression, and duplication were filtered out. The screening conditions for differential miRNAs were FDR < 0.05, logfc > | 2.5 |.

Target Gene Prediction of Differential miRNA

The target genes were predicted using the miWalk web tool, and the signal pathways of the identified genes were enriched and analyzed. The up-regulated and down-regulated genes were analyzed by KEGG and GO, respectively, using an R Language Clusterprofiler package.7

Construction of PPI Network

STRING () is a database for predicting protein–protein interaction (PPI) networks, which can predict downstream protein interactions. For the target genes of the selected differential miRNAs, protein interaction was predicted using STRING, and the position of the protein in the regulatory network was predicted by observing the number of nodes. The target gene network was constructed using Cytoscape software.

Survival Curve

The target genes were verified in the TCGA database, and the survival curve was drawn.

Statistical Methods

The differential genes were screened using an R Language Limma package, and the screening criteria were up-regulated or down-regulated with a difference multiple > 2.5 and P < 0.05. GO and KEGG analyses showed P < 0.05 to be statistically significant.

Results

Analysis of Differentially Expressed Genes

A total of 34 differentially expressed miRNAs were screened by comparing the NB group and the control group, of which 2 were up-regulated and 12 were down-regulated in the NB group (Table 1). The volcano map (Figure 1) and heat map (Figure 2) were drawn with differential genes. The color depth represents the level of gene expression, and the heat map can directly reflect the expression of differential miRNAs.

Table 1

Differentially Expressed miRNAs

	LogFC	P.Value	Change
hsa-miR-877-5p	6.035876	4.22E-05	UP
hsa-miR-7704	5.846523	0.000203	UP
hsa-miR-127-3p	3.32826	0.000308	UP
hsa-miR-369-5p	3.068657	0.00172	UP
hsa-miR-654-5p	3.769141	0.004878	UP
hsa-miR-1468-5p	3.233108	0.00564	UP
hsa-miR-132-5p	3.147397	0.006992	UP
hsa-miR-223-5p	3.394489	0.010416	UP
hsa-miR-412-5p	4.765532	0.011007	UP
hsa-miR-5010-5p	3.769372	0.018675	UP
hsa-miR-432-5p	3.430824	0.020249	UP
hsa-miR-1284	3.734122	0.020616	UP
hsa-miR-323a-3p	3.729744	0.021401	UP
hsa-miR-410-3p	3.78004	0.026546	UP
hsa-miR-382-5p	3.084026	0.029777	UP
hsa-miR-490-5p	3.598977	0.032409	UP
hsa-miR-487a-5p	3.327275	0.034281	UP
hsa-miR-744-3p	2.999005	0.035594	UP
hsa-miR-615-3p	2.690153	0.036106	UP
hsa-miR-4659a-3p	3.897584	0.038346	UP
hsa-miR-4435	3.353725	0.039823	UP
hsa-miR-3614-5p	3.56925	0.043763	UP
hsa-miR-1229-3p	−2.65241	0.000481	DOWN
hsa-miR-1292-5p	−3.61094	0.000721	DOWN
hsa-miR-6894-3p	−3.29068	0.001289	DOWN
hsa-miR-4685-3p	−4.07869	0.00293	DOWN
hsa-miR-200a-3p	−4.28626	0.019782	DOWN
hsa-miR-30a-5p	−3.45932	0.020238	DOWN
hsa-miR-671-5p	−2.51995	0.026266	DOWN
hsa-miR-1260b	−3.32059	0.034006	DOWN
hsa-miR-21-3p	−3.15339	0.038475	DOWN
hsa-miR-190b	−2.7182	0.040319	DOWN
hsa-miR-421	−2.53751	0.045866	DOWN
hsa-miR-29b-3p	−3.37678	0.046898	DOWN

Figure 1

Volcano map of differential miRNA.

Figure 2

Heat map of 22 upregulated miRNAs and 12 down-regulated miRNAs.

Prediction of Target Genes with Differentially Expressed miRNAs

Using the miWalk web tool, we predicted the target genes of 34 differentially selected miRNAs, including 1163 target genes predicted by 22 up-regulated miRNAs and 1474 target genes predicted by 12 down-regulated miRNAs.(See and for details).

KEGG Enrichment Analysis of Target Genes

The clustering profiler package of R language was used to analyze 1163 up-regulated target genes. The FOXO signaling pathway, mTOR signaling pathway, AMPK signaling pathway, and other pathways were observed in the results (Figure 3); a total of 1474 down-regulated target genes were analyzed by KEGG enrichment, and the results showed the relaxin signaling pathway, EGFR tyrosine kinase inhibitor resistance, and other pathways (Figure 4).

Figure 3

KEGG enrichment results of up-regulated target genes.

Figure 4

KEGG enrichment results of down-regulated target genes.

GO Enrichment Analysis of Target Genes

The cluster profiler package of R language was used for GO enrichment analysis of 1163 up-regulated target genes, and the results showed the physiological processes, such as axon formation, regulation of chemical synaptic transmitters, and activation of protein kinase activity, to be involved (Figure 5). The GO enrichment analysis of 1474 down-regulated target genes was performed, and the results showed such physiological processes as the regulation of synapses and the regulation of trans synaptic signals to be involved (Figure 6).

Figure 5

GO enrichment results of up-regulated target genes.

Figure 6

GO enrichment results of down-regulated target genes.

The target genes predicted by differential miRNAs were introduced into the STRING database, and the results were input to the Cyberscape software to construct the target gene network map. The top 16 genes were selected as the key genes (Figure 7).

Figure 7

16 key genes.

According to the data of NB in the TCGA plan, the key genes mentioned above were divided into a high-expression group and a low-expression group, and a survival curve was drawn. The results showed the survival curve of high- and low-expression groups of TP53 and ATM genes to be statistically significant (P < 0.05) (Figures 8 and 9).

Figure 8

ATM gene survival curve.

Figure 9

TP53 survival curve.

Discussion

NB is an embryonic tumor originating from neural crest stem cells, which can occur in the paravertebral sympathetic ganglion of adrenal medulla or neck, chest, abdomen, or pelvis. The most common clinical manifestations are painless abdominal masses.8 The diagnosis of NB should be combined with laboratory examination, radiological imaging, and pathology. The tumor would have often already metastasized by the time of diagnosis, and prognosis is often very poor. The high rate of metastasis, adverse prognosis, and unique clinical characteristics of advanced NB have inspired in-depth study of its biology and new methods of treatment.

miRNAs are small non-encoding RNA molecules about 22 nt in length. They are usually regulated by binding with their complementary mRNA.9 miRNAs participate in various aspects of the pathogenesis of NB, which strongly suggests that they have potential in the clinical management of NB. In this study, GSE128004 in the GEO database was analyzed using an R Language Limma package. The data from second-generation sequencing of the peripheral blood of 15 volunteers with NB and 3 healthy volunteers were obtained. A total of 34 miRNAs were obtained, including 22 up-regulated miRNAs and 12 down-regulated miRNAs. The target genes predicted by miRNAs were obtained by using the Miwok web tool, including 1163 up-regulated genes, 1474 down-regulated genes, and sub-genes. These genes were not enriched by KEGG or GO. Finally, 16 key genes were screened by PPI, including EGFR, TP53, KRAS, mapk1, AR, arb1, NEDD4L, asb11, Skp1, srsf1, HDAC1, ATM, PIK3CA, NR3C1, SRC, and snap25. TP53 is one of the most mutated genes in human malignant tumors, although it has a lower mutation rate in NB. Vid Mlakar et al10 showed that TP53 can induce the death of NB cells by regulating prima-1met through the glutathione level, and TP53 loses its antitumor function once it is inactivated.11,12 KRAS plays an important role in regulating cell proliferation, differentiation, and apoptosis. Domestic and international studies have shown that KRAS is highly expressed in NB;13,14 so, it can be used as the target of treatment. EGFR, PIK3CA, SRC, and other genes have been found to be highly expressed in NB,15–17 and this is consistent with this study. In addition, there are not many studies on AR, arb1, NEDD4L, and other genes, which can be used to direct future research.

The potential differences in miRNA levels in the peripheral blood of patients with NB and healthy volunteers were analyzed by bioinformatics. Out of the 34 miRNAs, 22 were up-regulated and 12 were down-regulated. The expression of Hsa Mir 490–5p and HSA mir-323a-3p in patients with NB was up-regulated, which was consistent with the findings reported by Zeka et al.18

KEGG analysis showed that the enrichment pathway of up-regulated genes mainly includes the FOXO signal pathway, mTOR signal pathway, AMPK signal pathway, and ErbB signal pathway. The enrichment pathway of down-regulated genes mainly includes tyrosine kinase inhibitor resistance, insulin secretion, and non-small cell lung cancer. FOXO signaling pathways control a variety of cell functions, such as cell death, metabolism, metastasis, and tumor angiogenesis. Stefan Saler et al19 confirmed that FoxO3 can promote angiogenesis and chemoresistance in vitro, and inhibition of FOXO signaling pathway can promote chemosensitivity in advanced NB,20 which may play an important role in the metastasis and treatment of NB in children. The mTOR signaling pathway is very important for cell growth, proliferation, and survival, and it is up-regulated in NB in children,21 which is consistent with this study. Some studies22 have shown that the ErbB signaling pathway plays an important role in various kinds of malignant tumors, including NB. Another study23 showed that signal transduction of ErbB receptor tyrosine kinase promotes NB in vitro and in vivo.

Through the analysis of GO enrichment, we found that the enrichment pathway of up-regulated genes mainly includes axon formation, chemical synaptic neurotransmitter modulation, and protein kinase activity. The enrichment pathway of down-regulated genes mainly includes cell morphology calibration, regulation of neuronal projection development, and gland development. During the differentiation of neural crest cells in the sympathetic adrenal system, any significant interruption or failure of regulatory function will lead to the transformation of neural crest cells into NB of the adrenal and sympathetic ganglia. This was mainly found to be related to the growth and development of neural cells in our GO analysis. Important processes included axogenesis, chemical synaptic neurotransmitter modulation, and regulation of neuronal projection development.

We also tested the key genes selected in the TCGA tumor database, and the results showed that the expression of TP53 and ATM genes is significant for survival. NB is the most common extracranial solid tumor in children. Its clinical manifestations are not specific. Metastasis has been found at the time of diagnosis, and prognosis is very poor. The 34 miRNAs and 16 key genes found by bioinformatics can play a role in the diagnosis of NB. These findings can provide some clues and ideas for molecular targeted treatment of NB and can predict the prognosis of the disease to a certain extent, but further experiments are still needed to verify these effects.