Joint analysis of the metabolomics and transcriptomics uncovers the dysregulated network and develops the diagnostic model of high-risk neuroblastoma

High-risk neuroblastoma (HR-NB) has a significantly lower survival rate compared to low- and intermediate-risk NB (LIR-NB) due to the lack of risk classification diagnostic models and effective therapeutic targets. The present study aims to characterize the differences between neuroblastomas with different risks through transcriptomic and metabolomic, and establish an early diagnostic model for risk classification of neuroblastoma.Plasma samples from 58 HR-NB and 38 LIR-NB patients were used for metabolomics analysis. Meanwhile, NB tissue samples from 32 HR-NB and 23 LIR-NB patients were used for transcriptomics analysis. In particular, integrative metabolomics and transcriptomic analysis was performed between HR-NB and LIR-NB. A total of 44 metabolites (P < 0.05 and fold change > 1.5) were altered, including 12 that increased and 32 that decreased in HR-NB. A total of 1,408 mRNAs (P < 0.05 and |log2(fold change)|> 1) showed significantly altered in HR-NB, of which 1,116 were upregulated and 292 were downregulated. Joint analysis of both omic data identified 4 aberrant pathways (P < 0.05 and impact ≥ 0.5) consisting of glycerolipid metabolism, retinol metabolism, arginine biosynthesis and linoleic acid metabolism. Importantly, a HR-NB risk classification diagnostic model was developed using plasma circulating-free S100A9, CDK2, and UNC5D, with an area under receiver operating characteristic curve of 0.837 where the sensitivity and specificity in the validation set were both 80.0%. This study presents a novel pioneering study demonstrating the metabolomics and transcriptomics profiles of HR-NB. The glycerolipid metabolism, retinol metabolism, arginine biosynthesis and linoleic acid metabolism were altered in HR-NB. The risk classification diagnostic model based on S100A9, CDK2, and UNC5D can be clinically used for HR-NB risk classification.

to January 2022.Inclusion criteria consisted of: (1) confirmed diagnosis of pathological NB; (2) the risk grade was clinically assessed based on COG classification; (3) informed consent was signed by children or their parents.Exclusion criteria consisted of: (1) complications of other diseases; (2) informed consent was not signed by children or their parents.Plasma samples were obtained from fasting plasma of NB patients on the morning of surgery and immediately frozen in a −80 °C refrigerator for metabolomics analysis.NB tissue samples were obtained from the patient's surgically resected tumor tissue and stored directly in liquid nitrogen for transcriptomic analysis.Tables S1 and S2 showed that there were no significant differences in age and gross tumor volume between HR-NB and LIR-NB whereas the tumor metastasis condition, MYCN amplification, and radiological risk factors were significantly different between HR-NB and LIR-NB samples.

Metabolomics analysis via high performance liquid chromatography-mass spectrometry (HPLC-MS)
The plasma samples were taken from −80 °C storage and immediately thawed in a 4 °C refrigerator.After 10 s of vortexing, 150 μL plasma was added to a 1.5 mL microcentrifuge tube and 450 μL acetonitrile was added at 4 ℃.After 5 min of vortexing at 3000 r/min followed by centrifugation at 13,000 r/min for 15 min (4 ℃), 300 μL supernatant was carefully extracted.Quality control (QC) samples were prepared by mixing equal amounts of supernatant from all samples and were used to evaluate the stability of the overall experimental results.All of the extracts were analyzed using an Agilent 6210 time-of-flight MS system with an Agilent 1100 HPLC, a photodiode array detector, a high-resolution-time-of-flight-MS with an electrospray ionization source and an Agilent workstation.The chromatographic separation was performed on an Agilent Poroshell 120 EC-C18 (2.7 μm, 3.0 × 100 mm) column.Metabolomics data were collected as follows: mobile phase: A = 0.1% formic acid water, B = 0.1% formic acid acetonitrile, elution conditions: 0-3 min, 5-60% B; 3-25 min, 60-90% B; 25-30 min, 90-100% B; 30-40 min, 100% B. Settings consisted of an injection volume of 10 μL, a column temperature of 30 °C and the flow rate of 0.3 mL/min.Mass spectrometry (MS) negative and positive mode conditions: nitrogen as drying gas, nitrogen temperature 325℃, flow rate 12 L/min, atomization pressure 35 psi; capillary voltage: positive mode 4,000 V, negative mode 3,500 V; fragmentation voltage: positive mode 215 V, negative mode 175 V, separator voltage 60 V; mass acquisition range: all the negative modes were 0.05-1.5 KDa.The samples were analyzed by HPLC-MS to obtain the original data files.Agilent Masshunter HPLC-MS software was used to convert the original data files into a common format.Based on R language platform, XCMS software package was used to conduct retention time (RT) calibration, peak recognition, noise filtering and peak matching for the obtained .mzDataformat files, and to set the allowable deviation of mass charge ratio and RT (mass/charge ratio tolerance = 0.025DA, RT tolerance = 0.5 min).The metabolites with RT deviation of 0.5 min and mass number deviation of 0.025 Da were considered to be the same metabolites.Finally, the data matrix containing mass/charge ratio, RT, peak area and other information were obtained.Metabolites were identified using both primary and secondary MS techniques.Firstly, the acquired primary MS information was subjected to targeted secondary MS analysis to obtain secondary MS information and provide a reference for subsequent qualitative analysis.Next, according to the precise mass number of excimer ions such as [M + H] + ions and the high-resolution target MS/MS spectrum, combined with the fragmentation laws of various metabolites, and through the online database METLIN (http:// metlin.scrip ps.edu/), HMDB (http:// hmdb.ca/) retrieval and literature retrieval methods, analysis and derivation of the possible structure of differential metabolites, get the information of candidate metabolites.The metabolomics analysis, partial least-squares discrimination analysis (PLS-DA), heatmap, volcano map, enrichment analysis, pathway analysis and biomarkers, were conducted by the MetaboAnalyst (https:// www.metab oanal yst.ca/ Metab oAnal yst/ home.xhtml).

Transcriptomics detection through RNA-sequencing (RNA-seq) analysis
Total RNA was extracted using TRIzol reagent according to the manufacturer's protocol.RNA purity and quantification were evaluated using the NanoDrop 2000 spectrophotometer (Thermo Scientific, USA).RNA integrity was assessed using the Agilent 2100 Bioanalyzer (Agilent Technologies, Santa Clara, CA, USA).The libraries were constructed using TruSeq Stranded mRNA LT Sample Prep Kit (Illumina, San Diego, CA, USA) according to the manufacturer's instructions.The transcriptome sequencing and analysis were conducted by OE Biotech Co., Ltd.(Shanghai, China).The libraries were sequenced on an Illumina HiSeq X Ten platform and 150 bp pairedend reads were generated.About 48.349 M raw reads for each sample were generated.Raw data (raw reads) of fastq format were firstly processed using Trimmomatic 18 and the low quality reads were removed to obtain clean reads.Around 47.459 M clean reads for each sample were retained for subsequent analyses.The clean reads were mapped to the human genome (GRCh38) using HISAT2 19 .Fragments per kilobase of exon model per million mapped fragments (FPKM) 20 of each gene was calculated using Cufflinks 21 and the read counts of each gene were obtained by HTSeq-count 22 .Differential expression analysis was performed using the DESeq (2012) R package 23 .P value < 0.05 and |log2(fold change)|> 1 were set as the threshold for significantly differential expression.Hierarchical cluster analysis of differentially expressed genes was performed to demonstrate the expression pattern of genes in different groups and samples.Open database sources, including the Gene Ontology (GO), Kyoto Encyclopedia of Genes and Genomes (KEGG) [24][25][26] , MetaboAnalyst, Human Metabolome Database and National Center for Biotechnology Information were used to identify metabolic pathways.

Joint analysis of the metabolomics and transcriptomics
Finally, comprehensive transcriptomics and metabolomics analyses were performed using the MetaboAnalyst 5.0.The joint-pathway analysis module was selected to perform topological analysis by entering official gene symbols and compound names with optional fold changes to assess the potential importance of individual molecules (i.e., nodes) according to their position in the network.

Detection the HR-NB plasma potential biomarkers using quantitative PCR based approach
The 24 differential genes were selected as potential candidate biomarkers based on transcriptomics results and literature.Corresponding primers for the differential genes were designed through the website of the National Center for Biotechnology Information (www.ncbi.nlm.nih.gov).Quantitative PCR (qPCR) was performed using the HiScript III All-in-one RT SuperMix kit and AceQ qPCR SYBR Green Master Mix kit (Vazyme, Nanjing) based on kit instructions.The housekeeping gene NAGK was selected as an internal control for mRNA abundance.Fold changes in the levels of target gene mRNA were determined using the formula 2 −ΔΔCt .

Develop the risk diagnostic model of NB using qPCR based approach
Logistic regression analysis was used to establish a regression equation for the test set (27 cases of HR-NB, 18 cases of LIR-NB) and validated against the validation set (15 cases of HR-NB, 15 cases of LIR-NB).SPSS software version 21.0 (IBM Corp., Armonk, New York) was used for data processing and GraphPad Prism 8.0 (GraphPad Software Inc., San Diego, USA) was used for mapping.

Results and discussion
The general idea of this research is shown in Fig. 1.Metabolomics analysis was performed on 96 plasma samples (58 cases of HR-NB, 38 cases of LIR-NB) and transcriptomics analysis was performed on 55 NB tissue samples (32 cases of HR-NB, 23 cases of LIR-NB).By integrating metabolomics and transcriptomics data, using PLS-DA, heatmap, enrichment analysis, pathway analysis and other methods for analysis and processing, the abnormal pathway network of HR-NB was comprehensively analyzed and potential clinical therapeutic targets were discovered.Finally, a diagnostic model was established, which is expected to be used for the early diagnosis of HR-NB.

The Metabolome differences between HR-NB and LIR-NB
To explore the differences in metabolites among HR-NB, IR-NB and LR-NB, plasma metabolomics analysis was first conducted using a non-targeted metabolomics-based approach.The PCA plot demonstrates the strong clustering of QC samples in both positive and negative modes, providing evidence for the robustness of our findings (Fig. S1).In order to comprehensively understand the metabolomics of NB with different risk levels, the PLS-DA among HR-NB, IR-NB and LR-NB was conducted.As shown in Fig. S2, the IR-NB group and the LR-NB group were clustered together and were situated far away from the HR-NB group, indicating that the LR-NB and IR-NB groups were closely related.However, the plasma metabolites in the IR-NB and LR-NB groups were significantly different from the HR-NB metabolites both in positive mode and negative mode, consistent with the clinical trend of high 5-year survival rates in children with LIR-NB [3][4][5] .Therefore, in order to better distinguish HR-NB, IR-NB, and LR-NB, IR-NB and LR-NB were merged into one group (LIR-NB).Following, the differences between LIR-NB and HR-NB metabolites were systematically studied.
In order to intuitively express metabolomic differences between LIR-NB and HR-NB, cluster analysis on the plasma metabolites of NB was conducted based on the correlation of compounds and presented in the form of a heatmap (Figure S3), illustrating the differences between the two groups.In order to understand the metabolomics of HR-NB and LIR-NB, the PLS-DA among HR-NB and LIR-NB was conducted first.In positive mode, cumulative R 2 Y was at 0.722 and Q 2 was at 0.575 (Fig. 2A) and in negative mode, cumulative R 2 Y was at 0.705 and Q 2 was at 0.506 (Fig. 2B), indicating that the model had excellent predictive ability.Volcano plots show the metabolites in the HR-NB and LIR-NB groups in positive and negative modes (Fig. 2C,D).Plasma metabolites with fold change > 1.5 and P < 0.05 in the volcano plot were identified as differential metabolites.Therefore, a total of 31 metabolites changed significantly in positive mode, including 21 up-regulated and 10 down-regulated   1).In negative mode, 13 differential compounds were identified, including 4 up-regulated compounds and 9 down-regulated compounds (Table 2).According to the differential metabolite results, the expression of phospholipids, amino acids and carnitine was significantly altered in HR-NB.Phospholipids play an important role in the composition and determination of biophysical properties of biological membranes, including their mobility, cell signaling, cell-cell interactions in tissues, and molecular transport 27,28 .Therefore, we believed that phospholipids played an important role in the development of HR-NB due to the marked alteration of phospholipids.Lysophosphatidylcholine (LPC) has been found to be a biomarker for several tumors, including prostate and lung cancer.LPC has also been associated with inflammation, oxidative stress, insulin resistance, apoptosis, lipid remodeling and signal transduction lipogenesis, which is consistent with our metabolomic findings of LPC alterations, thus we suggested that LPC played an important role in HR-NB [29][30][31][32][33] .Furthermore, in order to further visualize the differential metabolites, heatmaps of differential compounds were drawn according to the correlation of differential compounds (Fig. 2E,F).The heatmap shows the up-regulation of compounds such as LysoPc(0:0/18:0) and the down-regulation of metabolites such as SM(d18:2/14:0), indicating that these compounds have a certain correlation to HR-NB.Hence, a total of 44 differential metabolites were found in metabolomics, demonstrating the significant differences between HR-NB and LIR-NB.

The altered pathways and biomarkers between HR-NB and LIR-NB
In order to discover abnormal metabolic pathways based on the discovery of abnormal metabolites, enrichment analysis and pathway analysis were conducted.According to the 44 most altered metabolites, the mitochondrial β-oxidation of long chain saturated fatty acids, β-oxidation of very long chain fatty acids, betaine metabolism, carnitine synthesis, oxidation of branched chain fatty acids, plasmalogen synthesis, mitochondrial β-oxidation of short chain saturated fatty acids, methionine metabolism, fatty acid metabolism and glycine and serine metabolism were enriched in the enrichment analysis (Fig. 3A).Recent studies have identified fatty acid oxidation (FAO) as an important source of metabolites that promote cancer growth, showing high activity in various cancers such as breast cancer, glioma, ovarian cancer, which is consistent with the results of our enrichment analysis, FAO is also significantly altered in HR-NB, indicating the important role of FAO in HR-NB [34][35][36][37] .Interestingly, it has been shown that the activation of FAO can protect cancer cells from effects such as glucose deprivation and hypoxia whereas most cancer cells use glycolysis as the main energy source, therefore we deemed that the glycolysis in HR-NB activated FAO and affected the metabolic pathway of HR-NB.In order to further search for potential abnormal metabolic pathways and visualize the results, pathway analysis was performed.Mitochondrial β-oxidation of long-chain saturated fatty acids, β-oxidation of very long chain fatty acids, plasmalogen synthesis, betaine metabolism, oxidation of branched chain fatty acids and methionine metabolism were significantly changed in the pathway analysis (Fig. 3B).Therefore, 10 significantly altered metabolic pathways were identified by enrichment analysis and 6 significantly altered metabolic pathways were identified by pathway analysis, which are helpful to understand the abnormality of NB pathway network.
To explore the plasma biomarkers of HR-NB in metabolomics and provide a non-invasive strategy for NB risk classification, receiver operating characteristic (ROC) curves of differential metabolites were analyzed.The iconic biomarkers SM (d16:1/16:0) and SM (d18:2/14:0) were found (Fig. 3C,D).Other biomarkers are shown in Fig. S4, Table S3.Area under curve (AUC) of the ROC for all biomarkers were > 0.7, indicating that these metabolites may be potential biomarkers of HR-NB.Both SM (d16:1/16:0) and SM (d18:2/14:0) belong to the sphingolipids, a class of lipids that are expressed in all eukaryotic cells.They are relatively abundant in neuronal tissue cells, enriched in the plasma membrane of cells, and more strongly enriched in certain plasma membrane structural domains, such as the myelin sheath of oligodendrocytes or the parietal membrane of enterocytes 38 .In addition, they play a signaling role in many cellular functions, such as growth, cell cycle progression, differentiation and apoptosis 39,40 .Therefore, we inferred that alterations in sphingolipid content were associated with the development of NB and might be a valid biomarker for NB.Hence, 10 altered metabolic pathways were found in the enrichment analysis, 6 altered metabolic pathways were found in the pathway analysis, and 7 biomarkers were found by the biomarker analysis.These findings are helpful to understanding the abnormality of NB pathway network and targeted therapy.

Transcriptomics analysis uncovers the abnormal expression gene between HR-NB and LIR-NB
To further analyze the differences between HR-NB and LIR-NB, transcriptomics differences between 32 HR-NB tissues and 23 LIR-NB tissues were investigated.The detailed results of total RNA concentration, A 260 /A 280 , A 260 / A 230 , 28S/18S and RNA integrity number of the extracted samples are shown in Table S4.The RNA integrity number values extracted in this study were all greater than 7.The pre-processing results of sequencing data quality showed that the RawBases values of each sample ranged from 6.49 to 7.76 G, the CleanBases values ranged from 6.00 to 7.22 G, and the percentage of Q30 bases in each sample ranged from 92.59 to 95.18%.The GC content of each sample ranged from 44.29 to 50.53% (Table S5).Combining with FPKM (Figure S5) and the total number of mRNAs detected in the samples (Figure S6), the RNA quality of both groups met the standard and could be used for subsequent analysis.
In order to visually express the transcriptomics differences between LIR-NB and HR-NB, the tissue RNA of NB was clustered based on RNA correlation and presented it in the form of a heatmap (Fig. S7), illustrating the obvious differences between the two groups.NB tissue RNAs with P < 0.05 and |log 2 (fold change)|> 1 in the volcano plot (Fig. 4A) were identified as differential genes.A total of 1,408 differential genes were found, including 1,116 for the number of up-regulated genes and 292 for the number of down-regulated genes (Fig. 4B).The top 20 up-regulated and down-regulated genes in HR-NB and LIR-NB tissues are shown in Tables 3 and 4. It has been suggested that in FABP4-mediated macrophages can promote the proliferation and migratory phenotype of NB cells by inactivating the NF-κB-IL1α pathway through ubiquitination of ATPB, which increases the migration, invasion and tumor growth of NB cells 41 and in turn, is consistent with our results where the level of FABP4 increased 14.238-fold.Therefore, we surmised that highly expressed FABP4 may be associated with the risk of NB.In addition, the study by Zhu, Y et al. found that the expression level of UNC5D in the LR-NB group was significantly higher than that in the HR-NB group, suggesting that a high level of UNC5D expression was significantly associated with a good prognosis 42 , which was consistent with our results where UNC5D expression exhibited a significantly different fold change of 0.281.Therefore, we thought that low expression of UNC5D might be significantly associated with aggressive tumor behavior.The top 100 differential genes were displayed by a cluster analysis heatmap (Fig. 4C), which was used to more intuitively distinguish the differences between HR-NB and LIR-NB, indicating that there were significant differences between the two groups.RT-PCR is a standard method for measuring gene expression due to its high sensitivity and convenience for validating RNA-seq results.Eleven of our selected differential genes were validated by RT-PCR.The RT-PCR results were consistent with transcriptomic results, illustrating the reliability of our transcriptomics results (Fig. 4D).Therefore, 1,408 significantly different genes were found between HR-NB and LIR-NB by transcriptomics, illustrating the differences between the two groups and the results were validated by RT-PCR.

KEGG and GO analysis between HR-NB and LIR-NB in transcriptomics
To discover abnormal pathways on the basis of differential genes, GO analysis and KEGG analysis were performed.GO annotation analysis of the obtained differential genes was performed to analyze the metabolic pathways in which the genes are located to determine their possible biological functions.The GO enrichment analysis top 30 bar chart shown in Fig. 5A.GO classification bar chart is shown in the Figure S8A.Regarding biological processes, the top three regulated expression were extracellular matrix organization, chemokine-mediated signaling pathway, neutrophil chemotaxis.Regarding cellular component, the top three significantly regulated expressions were extracellular space, extracellular region, extracellular matrix.Regarding molecular function, the top three significantly up-regulated expressions were chemokine activity, extracellular matrix structural constituent, CCR chemokine receptor binding.In addition, KEGG analysis (Figs. 5B, S8B) was performed and cytokine receptor interaction was found to be the most altered pathway, with 69 mRNAs enriched.Cytokines act by binding to specific receptors in the plasma membrane of target cells.The exploration of cytokine receptor interaction is important for understanding the pathogenesis of various human diseases, especially cancers, as well as for identifying potential therapeutic targets.In this study, it was confirmed that cytokine receptor interaction plays an important role in HR-NB.Therefore, more research on cytokine receptor interaction provides a new perspective to understand the molecular mechanism of HR-NB.Relationships between differential genes were visualized by protein-protein interaction circle diagrams (Figure S9), illustrating the close connection between genes such as HBB, HBD, HBA1, HBG2, etc.Therefore, the complex pathways in the transcriptomics between HR-NB and LIR-NB were excavated by GO analysis and KEGG analysis, such as glycan biosynthesis and metabolism and lipid metabolism, which provided experimental support for the subsequent targeted therapy of NB.

Integrated transcriptomics and metabolomics analyses between HR-NB and LIR-NB
In order to facilitate the systematic study of NB by linking important metabolites and genes through shared metabolic pathways, joint-pathway analysis was used.Nine significantly altered pathways were revealed by integrated analysis of transcriptomics and metabolomics data (Table 5).These altered pathways can be visualized in Fig. 6A, which includes glycerolipid metabolism, retinol metabolism, arginine biosynthesis and linoleic acid metabolism with P-values < 0.05 and impact ≥ 0.5, indicating that dysregulation of these pathways may lead to increased COG grading in NB.The differentially expressed genes associated with these pathways are shown in Table 6.
According to the results of the joint analysis, three pathways were found to be of interest: glycine, serine and threonine metabolism; glycolysis or gluconeogenesis and glycerolipid metabolism, each of which with significant differences in gene expression.In the network diagram of glycine, serine and threonine metabolism (Fig. 6B), the expression of GXT, ALAS2, AOC3, BHMT, CHDH and SDS were up-regulated in HR-NB compared to LIR-NB, which was consistent with metabolomic results (decreased levels of glycine and threonine and marked  www.nature.com/scientificreports/alterations in amino acid metabolism), but inconsistent in terms of serine (increased levels of serine).Glycine, serine and threonine are important precursors for protein, nucleic acid and lipid synthesis and participate in carbohydrate metabolism pathways.Glycolysis provides 3-phosphoglycerate to promote serine biosynthesis.
Serine can be interconverted with glycine through serine hydroxymethyltransferase and excess serine biosynthesis drives tumorigenesis.Therefore we theorized that the activation of HR-NB glycolysis may promote serine synthesis whereas the interconversion of serine and threonine may cause glycine, serine, and threonine results to vary inconsistently [43][44][45][46][47] .In the glycerolipid metabolism pathway (Fig. 6C), the expressions of AGPAT2, DGAT2, DGKB, GPAM and PNPLA2 were affected.Glycerolipid, an abundant cellular lipid with physiological roles in energy metabolism and membrane structure, is synthesized by modifying the 3-phosphoglycerol backbone through acylation and dephosphorylation reactions [48][49][50] .GPAT, AGPAT, and DGAT are key genes in this pathway by regulating the glycerol-3-phosphate pathway to control intracellular glycerolipid levels enzymes.Studies have shown that overexpression of GPAT3 leads to increased formation of triacylglycerols, and AGPAT-2 is also overexpressed in various cancers where specific inhibition of AGPAT can induce apoptosis in cancer cells 51,52 .These studies are consistent with the increased expression of GPAT, AGPAT, and DGAT genes in joint pathway analysis.Therefore, changes in lipid metabolism are closely related to NB and glyceride metabolism may be an important pathway for the development of HR-NB.Finally, the expression of FBP1, HK2, PCK1 was up-regulated and the expression of PKLR was down-regulated in the altered glycolysis or gluconeogenesis pathway (Fig. 6D).
In contrast to normal cells relying on mitochondrial oxidative phosphorylation to release energy, most cancer cells use glycolysis as their primary energy source, which produces ATP faster than oxidative phosphorylation to power rapid cell division 53 .This is consistent with our results, indicating the abnormality of glycolysis in HR-NB.www.nature.com/scientificreports/Therefore, these abnormally expressed genes may be used as targets for clinical therapy of NB.Other altered pathways are shown in Figure S10.Nine abnormal pathways were found through combined metabolomics and transcriptomics analysis and three interesting pathways were further mined and speculated to play a role in expanding the scope of targeted therapy.

Classification of NB with selected transcriptome candidate biomarkers
Currently, the risk stratification methods for HR-NB are intricate, clinical presentations are atypical, imaging examinations pose diagnostic challenges, and histological diagnosis remains controversial.HR-NB necessitates comprehensive evaluation involving imaging studies, pathology assessments and other investigations, resulting in the majority of cases being detected at advanced stages.In order to improve the low rate of early diagnosis of HR-NB, a more efficient early diagnostic model was established as a supplement to existing early diagnostic methods.The 24 candidate genes screened according to the results of transcriptomics and literature were analyzed by RT-PCR in order to find biomarkers that could be used for diagnosis (Table S6).The results for individual candidate genes were calculated using Eq. 2 −ΔΔCt and the sensitivity and specificity were further calculated.However, the results showed that the areas under the ROC curve of the S100A9, CDK2 and UNC5D were 0.685, 0.564 and 0.809, respectively, indicating that the detection performance of individual biomarkers was limited (Figure S11).Therefore, a diagnostic model combining the three biomarkers was established and the regression equation of the diagnostic model Y = 1.495 − 0.510 X 1 (S100A9) − 0.713 X 2 (CDK2) + 0.647 X 3 (UNC5D) was obtained by logistic regression analysis.In the validation set, the sensitivity and specificity of the diagnostic model were 80% and 80%, respectively, the cutoff value was 0.025 (Fig. 7A) and the AUC of the ROC was 0.836, indicating the effectiveness of the diagnostic model (Fig. 7B).In addition, a close association between these 3 genes and NB was also found.Overexpression of S100A9 enhanced the proliferation, migration and invasion of NB cells [54][55][56][57] , suggesting that S100A9 may be involved in the development of NB tumors 58 .In this study, S100A9 increased 7.23-fold in RNA-seq, therefore the high expression of S100A9 may be related to the risk classification of NB.
The CDK2 gene encodes a member of the serine/threonine protein kinase family of proteins that can regulate cell cycle progression 59,60 , which is consistent with our results of altered glycine, serine and threonine metabolism.Zhu et al. found that the expression level of UNC5D mRNA in the LR-NB group was significantly higher than that in the HR-NB group, suggesting that high levels of UNC5D expression were significantly associated with a good prognosis 42 .Our results are consistent with their study, showing that UNC5D expression was decreased in the HR-NB group and increased in the LIR-NB group.It has also been indicated that low UNC5D expression was significantly associated with aggressive tumor behavior and overexpression of UNC5D significantly inhibited malignant cell behavior, including cell proliferation and migration as well as tumor growth 61 .Therefore, S100A9, CDK2 and UNC5D are viable biomarkers that can be utilized in combination as an early diagnostic model for predicting the plasma risk classification of NB.This approach holds promise for non-invasive, low-cost detection of NB at an early stage.

Conclusions
In this study, 96 clinical plasma samples and 55 clinical NB tissue samples were analyzed and 1,408 differential genes and 44 differential metabolites were identified.The metabolomics and transcriptomics characteristics of HR-NB patients were demonstrated and a comprehensive network analysis was carried out.Compared with LIR-NB, HR-NB showed significant differences in glycerolipid metabolism, retinol metabolism, arginine biosynthesis and linoleic acid metabolism.The combination of S100A9, CDK2 and UNC5D was subsequently selected as the risk stratification early diagnostic model for HR-NB.The area under the ROC was 0.837 and the sensitivity and specificity were both 80.0%, indiciated that the diagnostic model could be used for early diagnosis of HR-NB and can be therapeutic targets in the future.In summary, the dysregulated network was examined by joint analysis of metabolomics and transcriptomics and a diagnostic model for HR-NB was developed.

Figure 1 .Figure 2 .
Figure 1. Outline of research method.Metabolomics analysis was performed on plasma samples and transcriptomics analysis was performed on NB tissue samples.By integrating metabolomics and transcriptomics data, the abnormal pathway network of HR-NB was analyzed and potential clinical therapeutic targets were discovered.Finally, a diagnostic model was established.

Figure 5 .
Figure 5. GO and KEGG analysis between HR-NB and LIR-NB in transcriptomics.(A) GO analysis of biological processes, cellular components and molecular functions of up-and down-regulated genes.(B) Bubble chart of the KEGG pathways (top 20) that differentially expressed genes significantly involved in.The dot color represents the p-value and the dot size represents the number of differential genes.

Figure 6 .
Figure 6.Integrated transcriptomics and metabolomics analyses of NB metabolic pathways.(A) Joint-pathway analysis of differential genes and differential metabolism.(B) The glycine, serine and threonine metabolism pathway, (C) the glycerolipid metabolism pathway and (D) the glycolysis or gluconeogenesis pathway with altered significantly genes in HR-NB compared to LIR-NB.Significant overexpression in red, significant downexpression in green.

Figure 7 .
Figure 7. Establishment of HR-NB early diagnosis model.(A) The prediction accuracies by the S100A9, CDK2 and UNC5D in test phase and validation set are compared between HR-NB and LIR-NB.(B) AUC value for prediction of HR-NB based on 3 plasma genes (S100A9, CDK2, UNC5D), AUC = 0.836.

Table 1 .
Differential expressed metabolites in HR-NB vs. LIR-NB in positive mode.

Table 2 .
Differential expressed metabolites in HR-NB vs. LIR-NB in negative mode.NO MetabolitesMass-

Table 3 .
The top 20 genes significantly up-regulated in HR-NB vs. LIR-NB.

Table 4 .
The top 20 genes significantly down-regulated in HR-NB vs. LIR-NB.

Table 5 .
Joint analysis pathways of differential metabolites and genes.

Table 6 .
Related differentially expressed genes by joint-pathway analysis.