Microarray datasets identification of DEGs
To identify DEGs, we performed background correction and normalization processing on the NB microarray datasets GSE12460, GSE13136, GSE14880, and GSE16237. Filtered DEGs through the limma package in the R language, with P-value < 0.05 and | log2FC | > 1 as the filtering conditions. Among them, the GSE12460 obtained a total of 1479 DEGs, including 402 up-regulated and 1077 down-regulated DEGs, the GSE13136 obtained a total of 2778 DEGs, including 1223 up-regulated and 1555 down-regulated DEGs, the GSE14880 obtained a total of 2068 DEGs, including 535 up-regulated and 1533 down-regulated DEGs, GSE16237 obtained a total of 3842 DEGs, including 496 up-regulated and 3346 down-regulated DEGs. The volcanos were used to show the differential expression of genes in the microarray datasets (Fig 1A).
Identification of overlapping DEGs
To identify the overlapping DEGs, we used the tool Venny to analyze the four microarray datasets and took the overlapping part of the four microarray datasets. In this study, we identified a total of 431 overlapping DEGs, including 90 up-regulated DEGs and 341 down-regulated DEGs (Fig 1B).
GO and KEGG enrichment of overlapping DEGs
To further understand the functions and pathways involved in overlapping DEGs, we used the DAVID database for enrichment analysis. The results showed that in Biological Processes (BP), DEGs were significantly enriched in cell adhesion, peripheral nervous system development, extracellular matrix organization, muscle organ development, and cell proliferation, in Cellular Components (CC), DEGs Significantly enriched in the neuromuscular junction, neuron projection, cell surface, filopodium, and axon, regarding Molecular Function (MF), DEGs were significant in collagen binding, glycosaminoglycan binding, structural constituent of muscle, receptor binding, actin-binding Enrichment. In addition, KEGG analysis showed that DEGs were significantly enriched in Proteoglycans in cancer, cell adhesion molecules, MicroRNAs in cancer, Glycine, serine, and threonine metabolism, and Axon guidance (Fig 1C).
PPI network construction and module analysis
To establish the protein-protein interaction, we used the STRING database to construct a PPI network for overlapping DEGs, involving a total of 194 nodes and 401 edges (Fig 2A). The non-interacting proteins were filtered out by Cytoscape, and the algorithm analysis of MCODE was performed, and the module with the highest MCODE score = 9.0 was selected as the hub module (Fig 2B). Enrichment analysis through BiNGO plug-in showed that DEGs in the module were mainly enriched in ligase activity, zinc ion binding, acid-amino acid ligase activity, ligase activity, forming carbon-nitrogen bonds, transition metal-ion binding, male germline stem cell division, asymmetric cell division, leg morphogenesis, and ubiquitin ligase complex. The visualization results are shown in Fig 3A. In addition, KEGG enrichment analysis showed that these DEGs were related to Acute myeloid leukemia, Transcriptional misregulation in cancer, MicroRNAs in cancer, and Pathways in cancer. The visualization results are shown in Fig 3B.
miRNA prediction of target genes
To predict the miRNA related to the target genes, we used the miRNet database to perform predictive analysis, imported the prediction results into Cytoscape for visualization, and constructed a target gene-miRNA interaction network (Fig 4). The results showed that FBXO9 interacts with 42 miRNAs (for example, hsa-mir-329-3p), FBXO17 interacts with 40 miRNAs (for example, hsa-mir-124-3p), HECW2 interacts with 43 miRNAs (for example, hsa-mir-19a-3p), MIB2 interacts with 17 miRNAs (for example, hsa-mir-335-5p), RNF19B interacts with 117 miRNAs (for example, hsa-mir-17-5p), RNF213 interacts with 140 miRNAs (for example, hsa-mir-101-3p), TRIM36 interacts with 40 miRNAs (for example, hsa-mir-25-3p), TRIM71 interacts with 163 miRNAs (for example, hsa-let-7a- 5p), ZBTB16 interacts with 53 miRNAs (for example, hsa-mir-15a-5p). The miRNAs in the network were further ranked by degree method, and the top five miRNAs with the highest scores are shown in Table 2. In addition, we further validated the pathways in which these miRNAs were involved through the miRPathDB database, with P-value < 0.05 as the filtering conditions. The results showed that these miRNAs were significantly annotated to MicroRNAs in cancer, mTOR signaling, Proteoglycans in cancer, Pathways in cancer, and TGF-beta Signaling Pathway. Most of the pathways, including the five pathways listed earlier, were supported by experimental validation (especially strong experimental validation). The 20 significant pathways are shown in Table 3.
Transcription factor prediction of target genes
To predict the TFs related to the target genes, we used the NetworkAnalyst database to conduct predictive analysis and constructed the target gene-TF interaction network (Fig 5A). The results showed that FBXO9 interacts with 38 TFs (for example, KLF16), FBXO17 interacts with 13 TFs (for example, ZEB1), MIB2 interacts with 35 TFs (for example, SIX5), RNF19B interacts with 17 TFs (for example, EGR2), RNF213 interacts with 51 TFs (for example, IRF1), TRIM36 interacts with 10 TFs (for example, CBX8), TRIM71 interacts with 19 TFs (for example, SOX5), ZBTB16 interacts with 4 TFs (for example, WT1). Further GO enrichment analysis of TFs in the network showed that BP was mainly enriched in peptidyl-lysine modification, covalent chromatin modification, and histone modification, CC was significantly enriched in transcription regulator complex, RNA polymerase II transcription regulator complex, and transcription repressor complex. MF was mainly involved in DNA-binding transcription activator activity, DNA-binding transcription activator activity, RNA polymerase II-specific and DNA-binding transcription repressor activity, RNA polymerase II-specific. The visualization results are shown in Fig 5B. In addition, KEGG enrichment analysis showed that these TFs were related to Transcriptional misregulation in cancers, TGF-beta signaling pathway, Thyroid hormone signaling pathway. The visualization results are shown in Fig 5C.
Survival analysis of hub genes
To study the correlation between gene expression and the prognosis of NB patients, we further analyzed the data with a sample size of 649 through the R2 database. Finally, we screened out FBXO9, HECW2, MIB2, RNF19B, RNF213, TRIM36, and ZBTB16 among these eight DEGs. The genes closely related to the patient's prognosis were used as hub genes. Compared with NB patients with relatively low hub genes expression in tumor tissues, NB patients with relatively high hub genes expression have a higher overall survival rate (P < 0.001) (Fig 6). We further analyzed the correlation between hub genes expression and MYCN expression in NB tissues. The results showed that the expression of hub genes and MYCN in NB tumor tissues were negatively correlated (P < 0.001) (Fig 7A). In addition, the expression of hub genes in MYCN non-amplified tumor tissues was higher than that in MYCN-amplified tumor tissues (P < 0.001) (Fig 7B).
Diagnostic Effectiveness of Biomarkers for Neuroblastoma
To make the results more reliable, we used the microarray dataset GSE73517 to verify. Expression levels for MYCN, FBXO9, HECW2, MIB2, RNF19B, RNF213, TRIM36, and ZBTB16 are shown in the heatmap (Fig 8A). Consistent with the above results, FBXO9, HECW2, MIB2, RNF19B, RNF213, TRIM36, and ZBTB16 have significantly higher expression in MYCN-nonamplified NB cells (P < 0.01) (Fig 8B). We analyzed the correlation between seven effective biomarkers and MYCN. Correlation results are shown in Fig 8C. In addition, the microarray dataset GSE73517 was used to verify the diagnostic effectiveness of NB biomarkers through ROC analysis. AUC more than 0.800 was considered to have excellent specificity and sensitivity to diagnose NB. As shown in Fig 8D, the AUC values for FBXO9, HECW2, MIB2, RNF19B, RNF213, TRIM36, and ZBTB16 were 0.864 (95% CI, 0.796×0.932), 0.817 (95% CI, 0.739×0.898), 0.902 (95% CI, 0.832×0.973), 0.896 (95% CI, 0.835×0.958), 0.781 (95% CI, 0.677×0.884), 0.918 (95% CI, 0.866×0.970) and 0.881 (95% CI, 0.812×0.950).
Experimental verification of hub genes
RT-PCR technology was used to determine the mRNA levels of seven hub genes (FBXO9, HECW2, MIB2, RNF19B, RNF213, TRIM36, and ZBTB16) in MYCN non-amplified cell line AS and MYCN amplified cell line BE2. As shown in Fig 9, the seven identified hub genes were all highly expressed in MYCN non-amplified NB cell lines (P < 0.05), consistent with the bioinformatics analysis results.
MYCN regulates the expression of FBXO9, RNF19B, and TRIM36
To study whether MYCN regulates the hub genes, we designed MYCN siRNA and overexpression plasmids. MYCN siRNA was transfected into MYCN amplified BE2 cells, and MYCN overexpression plasmid was transfected into MYCN non-amplified AS cells. Then detected the changes in the expression of hub genes at the mRNA levels. The results showed that after MYCN siRNA was transfected into BE2 cells, the expression of FBXO9, RNF19B, and TRIM36 increased at the mRNA levels, which were consistent with the predicted results (Fig 10A). After the MYCN overexpression plasmid was transfected into AS cells, the expression of FBXO9, RNF19B, and TRIM36 decreased at the mRNA levels (Fig 10B). These data indicated that MYCN has a counter-regulatory effect on the expression of FBXO9, RNF19B, and TRIM36 in NB cells.