
==== Front
Discov Oncol
Discov Oncol
Discover Oncology
2730-6011
Springer US New York

39302522
1334
10.1007/s12672-024-01334-0
Analysis
Identification of prognostic biomarkers in neuroblastoma using WGCNA and multi-omics analysis
Ke Yuhan Kyh202333@163.com

12
Ge Wenliang gewl@ntu.edu.cn

12
1 grid.440642.0 0000 0004 0644 5481 Department of Pediatric Surgery, Affiliated Hospital of Nantong University, Nantong, 226001 China
2 https://ror.org/02afcvw97 grid.260483.b 0000 0000 9530 8833 Department of Pediatric Surgery, Medical School of Nantong University, Nantong, 226001 China
20 9 2024
20 9 2024
12 2024
15 4698 7 2024
10 9 2024
© The Author(s) 2024
2024
https://creativecommons.org/licenses/by-nc-nd/4.0/ Open Access This article is licensed under a Creative Commons Attribution-NonCommercial-NoDerivatives 4.0 International License, which permits any non-commercial use, sharing, distribution and reproduction in any medium or format, as long as you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons licence, and indicate if you modified the licensed material. You do not have permission under this licence to share adapted material derived from this article or parts of it. The images or other third party material in this article are included in the article’s Creative Commons licence, unless indicated otherwise in a credit line to the material. If material is not included in the article’s Creative Commons licence and your intended use is not permitted by statutory regulation or exceeds the permitted use, you will need to obtain permission directly from the copyright holder. To view a copy of this licence, visit http://creativecommons.org/licenses/by-nc-nd/4.0/.
Background

Neuroblastoma (NB) is one of the most frequent parenchymal tumors among children, with a high degree of heterogeneity and wide variation in clinical presentation. Despite significant therapeutic advances in recent years, long-term survival in high-risk patients remains low, emphasizing the urgent need to find new biomarkers and construct reliable prognostic models.

Methods

In this study, data from neuroblastoma samples in the ArrayExpress database were utilized to identify key gene modules and pivotal genes associated with NB prognosis by weighted gene co-expression network analysis (WGCNA). The Gene Ontology (GO) and Kyoto Encyclopedia of Genes and Genomes (KEGG) functional enrichment analysis was performed using the DAVID database. Based on these hub genes, survival prognosis models were constructed and validated on an independent validation set in the Gene Expression Omnibus (GEO) database. Differences in biological functions and immune microenvironments and the sensitivity to pharmacological and immunotherapeutic treatments of patients in the high- and low-risk groups were examined by gene set enrichment analysis (GSEA) and immune infiltration analysis.

Results

WGCNA identified 14 gene modules and screened the module with the highest relevance to the International Neuroblastoma Staging System (INSS), containing 60 pivotal genes. GO and KEGG analyses demonstrated that these pivotal genes were mainly implicated in biological processes and signaling pathways including DNA replication, cell division, mitotic cell cycle, and cell cycle. Based on Lasso regression and COX regression analysis, a prognostic model containing DHFR, GMPS and E2F3 was constructed, and the RiskScore was significantly correlated with the 1-, 3- and 5-year survival of the patients. GSEA and immune infiltration analyses revealed significant differences in the levels of cell cycle-related pathways and immune cell infiltration between the high and low RiskScore groups. In particular, patients in the high-risk group are less likely to benefit from immunotherapy and may be better suited for treatment with drugs such as Oxaliplatin and Alpelisib.

Conclusion

This research systematically identified biomarkers related to NB prognosis and developed a reliable prognostic model applying WGCNA and multiple bioinformatics methods. The model has important application value in predicting patients’ prognosis, evaluating drug sensitivity and immunotherapy effect, and provides new ideas and directions for precise treatment of neuroblastoma.

Keywords

Neuroblastoma
Prognostic model
WGCNA
Biomarkers
Drug sensitivity
Precision medicine
Jiangsu Provincial Research HospitalYJXYY202204 Ge Wenliang issue-copyright-statement© Springer Science+Business Media, LLC 2024
==== Body
pmcIntroduction

Neuroblastoma (NB) is a type of embryonic tumor associated with the autonomic nervous system, suggesting that its origin lies in a precursor cell that is still in development and not yet fully specialized, originating from neural-crest tissues [1]. As one might anticipate with a condition affecting developing tissues, cases of neuroblastoma typically arise in infants and toddlers; the average age at the time of diagnosis is 17 months [2, 3]. NB being highly heterogeneous, its clinical presentation ranges from spontaneous regressions to highly invasive metastatic lesions resulting in a wide variation in prognosis [4]. Despite significant advances in the treatment of NB in recent years, particularly in high-risk patients where the 5-year survival rate has increased from less than 20% to more than 50%, the overall prognosis remains unsatisfactory [5]. Long-term survival in high-risk NB patients remains low, and treatment-related late effects and chronic health problems significantly affect the quality of life of survivors. Current treatments include surgery, chemotherapy, and radiotherapy, but there are no effective curative therapies for patients with recurrence, emphasizing that significant clinical challenges remain [4]. Previous study showed that despite improved survival in patients with high-risk neuroblastoma, the majority of patients with recurrence still do not have curative therapies and survivors face treatment-related late effects in the long term [6]. Krystal and Foster emphasized that, through decades of clinical trials and research, advances in treatments including surgery, cytotoxic chemotherapy, radiotherapy, and autologous stem cell transplant and radiotherapy have significantly improved the prognosis of patients with NB [7]. However, these treatments have had limited success in improving long-term survival.

Biomarker discovery based on sequencing data has become an important component in precision medicine [8]. In cancer research, biomarkers are used not only for diagnosis but also for predicting treatment response and monitoring disease progression. For example, Zhu et al. noted that biomarkers identified from large-scale transcriptomic data can be used for pan-cancer diagnosis and prognosis, facilitating the application of precision medicine for cancer treatment [9]. Sarhadi and Armengol noted that biomarkers discovered through sequencing data can be used in precision medicine for cancers including risk assessment, diagnosis, prognostic evaluation as well as prediction and monitoring of treatment response [10]. Shaikh et al. predicted potential biomarkers by computational methods such as text mining, network analysis, next-generation technologies, machine learning, and deep learning in a breast cancer study [11]. Demir Karaman and Işık identified prognostic biomarkers for multiple cancer types by analyzing multi-omic data markers [12]. Addanki et al. noted that liquid biopsies, including circulating tumor DNA (ctDNA) and circulating tumor cells (CTCs), are biomarkers identified from sequencing data that can be used for precision medicine in cancer [13].

This study analyzed neuroblastoma sample data from the ArrayExpress database using WGCNA and identified key gene modules and hub genes related to NB prognosis. Based on these hub genes, a survival prognosis model was constructed. Through these analytical methods, biomarkers related to NB prognosis were systematically identified, providing new ideas and directions for subsequent clinical applications and the improvement of personalized treatment plans.

Materials and methods

Data sources

NB data with number E-MTAB-8248 were collected from ArrayExpress database (https://www.ebi.ac.uk/biostudies/arrayexpress) as a training set, which included a total of 223 cancer samples. Meanwhile, the neuroblastoma dataset with the number GSE49710 was retrieved from Gene Expression Omnibus (GEO, https://www.ncbi.nlm.nih.gov/geo/) database as a test set, containing a total of 498 cancer samples.

Weighted gene co-expression network analysis (WGCNA)

WGCNA was conducted using the WGCNA package in R [14]. The INSS stage was used as a trait to perform WGCNA analysis on the training set, and the soft threshold intercept height was set to 0.9 to obtain the optimal soft threshold β = 7. The minimum number of genes was 100 in a module and the module merge height was 0.2. The module with the highest correlation was defined as the specificity module.

Functional enrichment analysis

The hub genes identified from the WGCNA obtained from the screening were analyzed for functional enrichment by the Gene Ontology (GO) and Kyoto Encyclopedia of Genes and Genomes (KEGG) using the Database for Annotation, Visualization, and Integrated Discovery (DAVID, https://david.ncifcrf.gov/) [15, 16]. Functional entries and pathways that were significantly enriched (Count ≥ 2, P Value < 0.05) were screened.

Construction of prognostic model

Univariate COX regression analysis enables the variables in the model to be closely associated with the clinical outcomes of the patients, thus improving the clinical relevance of our constructed model. Thus, the hub genes obtained from the WGCNA screening of the combined survival data were analyzed by univariate COX regression to screen the genes significantly associated with survival. Subsequently, we used Lasso regression to further compress the screened genes to remove redundant and irrelevant variables, thereby reducing the risk of overfitting the model. Lasso regression was performed using the R package glmnet [17], and the best model was determined by tenfold cross-validation. Finally, 3 genes associated with NB prognosis were identified by stepwise multifactorial regression analysis: DHFR, GMPS, and E2F3. The risk score formula was RiskScore = 0.522 * DHFR + 0.569 * GMPS + 0.408 * E2F3.

Model validation

Each sample in the training set was assigned with a risk score and receiver operating curve (ROC) analysis was conducted using the R package pROC to analyze the classification efficiency of prognostic prediction at 1, 3 and 5 years, respectively [18]. The samples were categorized into high-risk and low-risk groups, the overall survival of the two groups was compared, and the robustness of the model was verified in the GEO test set.

Gene set enrichment analysis (GSEA)

GSEA (FDR < 0.05) for the training set samples was conducted applying GSEA_4.2.2 software. The c2.cp.kegg (http://www.gsea-msigdb.org/gsea/msigdb/download_file.jsp?filePath=/resources/msigdb/7.4/c2.cp.kegg.v7.4.symbols.gmt) and c5.go gene sets (http://www.gsea-msigdb.org/gsea/msigdb/download_file.jsp?filePath=/resources/msigdb/7.4/c5.go.v7.4.symbols.gmt) were obtained from the Molecular Signatures Database portal (MSigDB, https://www.gsea-msigdb.org/gsea/msigdb/), and KEGG and GO gene set enrichment analyses were performed to analyze the differences in relevant biological functions and signaling pathways between the two risk groups.

Immune infiltration analysis

CIBERSORT is an algorithm that can be used to infer immune cell composition from complex tumor samples, which analyzes the proportions of different immune cells in the tumor microenvironment to understand the role of the immune system in tumorigenesis and progression [19]. The ESTIMATE algorithm is based on the “estimate” package and is used to estimate the stromal and immune cell scores of a tumor tissue sample, which in turn reflects the relative amounts of immune and stromal components in the sample [20]. We utilized both of these algorithms in order to assess differences in the extent of immune cell infiltration in patients from different NB risk groups. In addition, a spearman correlation-based analysis was performed in order to analyze the correlation between immune infiltration scores and prognostic biomarkers and risk scores.

Drug sensitivity and immunotherapy correlation

Drug sensitivities of patients in the training set were predicted using the R package oncoPredict [21], and the half-maximal inhibitory concentration (IC50) of drugs was calculated to screen for drugs significantly related to the risk score. Differences in immunotherapy efficacy between the two risk groups were assessed using TIDE software (http://tide.dfci.harvard.edu/) [22].

Results

WGCNA results

Firstly, Weighted gene correlation network analysis was performed in E-MTAB-8248 cohort.In constructing the scale-free network, the soft threshold intercept height was set to 0.9, at which time the optimal soft threshold β = 7 (Fig. 1A). After integrating similar gene modules, 14 gene modules were obtained (Fig. 1B). INSS Stage was analyzed as the analyzed trait with the main trait ME of each gene module for pearson correlation analysis. MEturquoise showed the highest and most significant positive correlation with INSS Stage (r = 0.32, p = 9e–07) with 3621 genes within the module (Fig. 1C). With gene significance (GS) ≥ 0.2 and |module membership (MM)|≥ 0.8 as screening conditions, 60 significant hub genes were finally obtained (Fig. 1D).Fig. 1 Results of Weighted gene correlation network analysis. A Soft threshold screening graph, screening the minimum value of 7 with R^2 higher than 0.9 as the soft threshold for constructing the topology network. B Gene-module control plot, the upper half of the plot is a gene tree, each branch represents a gene, and each gene corresponds to a module before and after merging below. C Module-trait correlation plot, the horizontal coordinates of the plot are shapes, the vertical coordinates are modules, the numbers in the squares are uncorrelation coefficients, and the numbers in the parentheses are the significance P-values, with the red color representing a positive correlation and the blue color representing a negative correlation. D GS-MM plot for green module, horizontal coordinate MM (module membership), the correlation of each gene with the trait gene; vertical coordinate GS (gene significance), the correlation of genes and traits within the module

Biological functions of hub genes

Towards the comprehension of Hub gene correlative biological functions and pathways, the obtained 60 genes were processed for GO and KEGG functional enrichment analysis, and screened for significantly enriched functional terms and pathways (Count ≥ 2, P value < 0.05). Hub genes were largely implicated in the GO Biological Process (BP) terms of mitotic spindle assembly checkpoint, cell division, mitotic spindle organization, mitotic cell cycle, and cell cycle, which are GO Biological Process terms (Fig. 2A). Figure 2B shows that the cellular components (CC) enriched for the central genes we screened include kinetochore, nucleoplasm, and nucleus, among others. In addition, molecular function (MF) enrichment analysis revealed that these genes were also enriched in multiple pathways such as ATP binding, DNA binding, and microtubule binding (Fig. 2C). Finally, KEGG enrichment analysis also demonstrated that the hub genes we screened were also significantly enriched in the cell cycle pathway (Fig. 2D). These results further reveal that these pathways have an impact on both the onset and development of NB.Fig. 2 Biological functions of hub genes. A Bubble diagram of the top 10 terms of biological process. B Bubble chart of the top 10 terms of cellular components. C Bubble chart of top 10 terms of molecular function. D Bubble plot of KEGG significance in top 10 terms. The horizontal coordinate in all plots is GeneRatio, which is the ratio of genes enriched to the total number of genes entered in the entry, the size of the circle in the plot represents the number of genes enriched to the entry, and the color represents the significance of the terms, with blue to red representing an increase in significance

Biomarker-based prognostic modeling to predict NB prognosis

In order to narrow down the scope of the model, the 60 genes obtained from WGCNA combined with the survival data from the E-MTAB-8248 cohort were subjected to univariate COX analysis, which showed that all 60 genes were closely related to survival. These 60 prognostically significant genes were further compressed using Lasso regression to shrink the gene number in the prognostic model. After tenfold cross-validation for model construction, the model was optimized when lambda = 0. 037, at which time 6 genes were included in the model (Fig. 3A). Finally, three prognosis-related genes in NB, DHFR, GMPS, and E2F3, were identified after multivariate COX analysis. The model was formulated as RiskScore = 0.522*DHFR + 0.569*GMPS + 0.408*E2F3. In the E-MTAB-8248 cohort, RiskScore predicted that the 1-year, 3-year, 5-year prognosis of patients had AUC values greater than 0.8 (0.89, 0.82, and 0.82, respectively). In the E-MTAB-8248 cohort, patients were categorized into a high RiskScore group (102 samples) and a low RiskScore group (121 samples) based on the optimal grouping intercept method of the survminer package. The survival rate and time of the high RiskScore patients were more unfavorable than the low RiskScore patients (Fig. 3C, D). The prognostic value of RiskScore was validated in validation set GSE49710 grouped in the same way. Evidently, the results in the validation set GSE49710 followed the trend of the results in the training set.The AUC values of RiskScore for predicting patients' 1-year, 3-year, and 5-year prognosis in GSE49710 were 0.86, 0.85, and 0.83, respectively, which were also all greater than 0.8, demonstrating that the RiskScore presented excellent predictive value (Fig. 3E). Similarly, patients in the high RiskScore had overall more unfavorable survival rates and times when compared to the low RiskScore group (Fig. 3F, D).Fig. 3 Biomarker-based prognostic model predicting the prognosis of NB. A LASSO COX analysis. The left panel shows a plot of the Lasso penalty term parameters, with log(lambda) values in the horizontal coordinates and degrees of freedom in the vertical coordinates, representing the cross-validation error, and the location where the cross-validation error is expected to be minimized in the actual analysis. The right figure shows the Lasso regression coefficients plot, with log(lambda) in the horizontal coordinate and the coefficients of the genes in the vertical coordinate, demonstrating the change of the coefficients of the different variables with the λ-penalty, when min Lambda = 0.037, according to which the position (min Lambda) determines the topmost cross-coordinate of the log (Lambda), and the top shows the number of featured genes6. B ROC curve of RiskScore predicting the prognosis of NB patients in the training set (E-MTAB-8248) model. C KM survival curves of patients in the high and low RiskScore groups in the training set (E-MTAB-8248) model. D Scatterplot of the distribution of RiskScore values and survival status of NB samples in the training set (E-MTAB-8248). E: ROC curve of RiskScore predicting the prognosis of NB patients in the model of the training set (GSE49710). F KM survival curves of patients in the high and low RiskScore groups in the training set (GSE49710) model. G Scatterplot of the distribution of RiskScore values and survival status of NB samples in the training set (GSE49710)

RiskScore was an independent prognostic factor for NB

Age and RiskScore were significantly correlated with OS of NB and acted as risk factors for prognosis of NB, as shown by the results of univariate COX analysis (Fig. 4A). Multivariate COX analysis showed that RiskScore could serve as an independent factor for NB prognosis (Fig. 4B). In addition, we observed the distribution of patients’ age, stage(inss), and tert_expression in high and low RiskScore subgroups. Evidently, patients with NB in the high RiskScore subgroup had higher age, higher stage (inss), and exhibited higher tert expression (Fig. 4C–E).Fig. 4 RiskScore is an independent prognostic factor for NB. A Forest plot of univariate COX analysis combining age, stage(inss), tert_expression, and riskscore. B Forest plot of multivariate COX analysis combining age, riskscore. C Boxplot of age information for NB patients in high and low RiskScore subgroups. D Boxplot of age information for NB patients in high and low RiskScore groupings. E Box line plot of age information for NB patients in high and low RiskScore groupings

Abnormal signaling pathways in high and low RiskScore groups

Abnormal signaling pathways in high and low RiskScore groupings were compared using the c2.cp.kegg gene set and the c5.go genome as benchmarks. It could be visualized that DNA REPLICATION, CHROMOSOME ORGANIZATION, REGULATION OF DNA REPLICATION, DNA REPAIR, RIBOSOME BIOGENESIS, SPLICEOSOME in the high RiskScore grouping, NUCLEOTIDE EXCISION REPAIR, CELL CYCLE, BASE EXCISION REPAIR, and HOMOLOGOUS RECOMBINATION were more active (Fig. 5). These pathways are cell cycle-related pathways, suggesting that activation of cell cycle pathways in patients with high RiskScore results in accelerated rates of cell division and proliferation.Fig. 5 GSEA results for samples in high and low RiskScore groupings

Differences in immune infiltration in RiskScore groupings

There were significant differences in the infiltration of five immune cells in different RiskScore groupings, Monocytes, T cells CD4 memory resting, Dendritic cells resting, T cells gamma delta, and Neutrophils (Fig. 6A). RiskScore showed a significant positive correlation trend with Plasma cells and Neutrophils, and RiskScore was significantly negatively correlated with T cells gamma delta and T cells CD4 memory resting (Fig. 6B). Meanwhile, DHFR, GMPS, and E2F3 showed consistency in correlating with the level of 22 immune cell infiltration, mainly showing positive correlation trend with Neutrophils, and negative relation to T cells CD4 memory resting, and T cells gamma delta (Fig. 6B). We also observed that the low RiskScore subgroup exhibited higher ImmuneScore and ESTIMATEScore overall (Fig. 6C). Finally, we observed that the expression levels of RiskScore, DHFR, GMPS, and E2F3 all showed a significant negative correlation with ImmuneScore (Fig. 6D).Fig. 6 Differences in immune infiltration in RiskScore subgroups. A Box line plot of immune cell infiltration of 22 species in RiskScore subgroups. B Heatmap of the correlation between the expression levels of RiskScore, DHFR, GMPS, and E2F3 and the level of immune cell infiltration of 22 species. C StromalScore, ImmuneScore, ESTIMATEScore of samples in the RiskScore grouping. D Heatmap of correlation between RiskScore, DHFR, GMPS, and E2F3 expression levels and StromalScore, ImmuneScore, and ESTIMATEScore

Assessment of chemotherapy and immunotherapy indications for patients in RiskScore subgroups

We also discussed the therapeutic sensitivity of chemotherapeutic small molecule drugs in high and low RiskScore subgroups. We found a significant correlation between 133 drugs and riskScore (|cor|≥ 0.3, p < 0.05). We show the top 10 drugs with positive and negative correlations, with major positive correlations including Doramapimod_1042, BMS-754807_2171, SB216763_1025, KU-55933_1030, NU7441_1038, AZD8055_1059, AZD6482_2169, JQ1_2172, Axitinib_1021, PF-4708671_1129, and major negatively related drugs include Dihydrorotenone_1827, Gallibiscoquinazole_1830, Fulvestrant_1816, Acetalax_1804, GDC0810_1925, BDP-00009066_1866, Alpelisib_1560, OF-1_1853, Oxaliplatin_1089, VE-822_1613 (Fig. 7A). Specifically, patients with low RiskScore were Axitinib_1021, SB216763_102, SB216763_1025, NU7441_1038, Doramapimod_1042, PF-4708671_1129, AZD6482_2169, and JQ1_2172 drug therapy More sensitive. Patients with high RiskScore had higher sensitivities to Oxaliplatin_1089, Alpelisib_1560, VE-822_1613, Acetalax_1804, Fulvestrant_1816, Gallibiscoquinazole_1830, OF-1_1853, GDC0810_1925 drug treatments were more sensitive (Fig. 7B). In addition, we noted that patients in the high RiskScore group exhibited higher TIDE scores, suggesting that patients therein may not be suitable for immunotherapy (Fig. 7C).Fig. 7 Assessment of chemotherapy and immunotherapy indications in RiskScore grouped patients. A Lollipop chart showing the 20 drugs most significantly correlated with riskscore, with drugs on the left side of the horizontal line showing a significant negative correlation with riskscore and drugs on the right side of the horizontal line showing a significant positive correlation with riskscore. B IC50 values of the 20 drugs in the high and low riskscore groups. C Box line plot of TODE scores of patients in the riskscore group

Discussion

This study successfully classified key gene modules related to NB by WGCNA and further analyzed the role of these genes in disease progression. Particularly, the MEturquoise module, whose high correlation with Stage (INSS) (r = 0.32, p = 9e–07) showed the importance of these genes in the development of NB. Subsequently, by systematically identifying key molecular markers related to NB prognosis, reliable prognostic models will be constructed, which will in turn better promote the precision treatment of NB.

In NB research, the identification and functional resolution of hub genes is an important step in understanding tumor mechanisms. The results of GO and KEGG functional enrichment analyses showed that hub genes were mainly involved in several key biological processes and signaling pathways. Cell division is the fundamental process of tumor cell proliferation. Tumor cells form masses and invade surrounding tissues by continuously undergoing uncontrolled cell division. The CDC20 gene promotes prostate cancer cell growth and invasion through stabilizing β-catenin in the cancer stem cell-like cells [23]. Mitotic Spindle Organization in mitosis plays a key role in ensuring proper chromosome segregation. Abnormalities in spindle organization may lead to chromosomal instability, which is a critical mechanism in cancer development. Studies have shown that the mechanism of spindle multicenter aggregation in cancer cells is closely correlated with the three-dimensional structure and invasive ability of cancer cells [24]. Uncontrolled cell cycle regulation is one of the main causes of tumorigenesis. Cancer cells evade normal cell growth control mechanisms by altering the function of cell cycle checkpoints, resulting in uncontrolled proliferation.NuMA protein establishes spindle bipolarity during mitosis by organizing microtubule astrosomes, which is essential for the normal progression of cell division [25]. KEGG analysis showed that the Hub genes were primarily involved in cell cycle, DNA replication and cancer signaling pathways. These pathways play important roles in tumorigenesis and progression. Abnormalities in cell cycle and DNA replication pathways can lead to unlimited cell proliferation, whereas activation of cancer signaling pathways promotes tumor cell survival, proliferation, and invasion. Overexpressed Cyclin D1 gene plays a critical partin the adaptation of cancer cells to the mechanical stresses within the tumor [26]. Taken together, Hub genes in neuroblastoma promote tumorigenesis and progression by regulating several key biological processes and signaling pathways. These findings could deepen the current understanding of neuroblastoma mechanism, providing potential targets for developing new strategies for NB treatment.

In prognostic models of NB, the DHFR, GMPS, and E2F3 genes play key roles. DHFR is involved in DNA synthesis and repair by catalyzing the reduction of dihydrofolate to tetrahydrofolate in cellular proliferation, and its high expression has been associated with many types of cancers, including breast and lung cancer. For example, DHFR is a prognostic factor dependent in ovarian cancer, and low levels of DHFR expression were associated with longer patient survival. In addition, inhibition of DHFR modulates anti-tumor immunity, showing a potential of using DHFR inhibitors in combating immunotherapeutic resistance in non-small cell lung cancer [27]. The role of GMPS in guanosine monophosphate synthesis and its aberrant expression are closely related to tumor growth. It has been shown that GMPS is part of a prognostic model constructed based on metabolism-related genes in pancreatic cancer [28]. In ovarian cancer, GMPS, PR, CD40 and p21 were identified as important prognostic markers [29]. E2F3 is a transcription factor that promotes rapid tumor cell growth by regulating cell cycle progression and DNA synthesis. E2F3 is overexpressed in many types of cancers, including prostate and bladder cancers.Higher expression of E2F3 in colon cancer was significantly associated with worse overall survival, suggesting that it may serve as a prognostic biomarker [30]. In addition, E2F3 acts as a promoter in colon cancer progression by activating the STAT3 pathway and accelerating the stemness of colon cancer cells [31]. Therefore, these three genes identified in our study may be important driving causative agents of NB, and their high expression may be associated with poor patient prognosis and serve as potential therapeutic targets. The RiskScore formed by the combination of these genes is able to predict the 1-, 3-, and 5-year survival rates of patients, which has important clinical applications. Prognostic models constructed using these genes can more accurately evaluate patients’ prognosis and help physicians better develop personalized treatment plans.

Conclusion

In conclusion, the roles of DHFR, GMPS, and E2F3 in neuroblastoma and their applications in prognostic models provide new perspectives and methods for precision tumor therapy, promoting a better understanding of the molecular mechanisms of NB, facilitating the development of new therapeutic strategies.

Limitations

The present study still has limitations, although the hub genes in Neuroblastoma were identified and a prognostic model was constructed by the WGCNA method, this is a predictive study based on computerized methods with sequencing data from patient samples. The performance of the Neuroblastoma prognostic model still needs to be validated with large samples of data in multiple clinical centers. Clinical value of the model. Secondly, the hub genes identified in this study have been less studied in Neuroblastoma, and their functions are still unclear, and their functions in Neuroblastoma still need to be explored through cellular and animal experiments.

Abbreviations

NB Neuroblastoma

WGCNA Weighted gene co-expression network analysis

GSEA Gene set enrichment analysis

CTCs Circulating tumor cells

CtDNA Circulating tumor DNA

GEO Gene Expression Omnibus

GO Gene Ontology

KEGG Kyoto Encyclopedia of Genes and Genomes

ROC Receiver operating curve

MSigDB Molecular Signatures Database

IC50 Half-maximal inhibitory concentration

Acknowledgements

Not applicable.

Author contributions

All authors contributed to this present work: [YHK] and [WLG] designed the study, [YHK] and [WLG] acquired the data. [YHK] and [WLG] improved the figure quality. [YHK] and [WLG] drafted the manuscript, [YHK] and [WLG] revised the manuscript. All authors read and approved the manuscript.

Funding

The study was supported by Jiangsu Provincial Research Hospital (YJXYY202204).

Data availability

The datasets generated and/or analyzed during the current study are available in the [GSE49710] repository, [https://www.ncbi.nlm.nih.gov/geo/query/acc.cgi?acc = GSE49710].

Declarations

Ethics approval and consent to participate

Not applicable.

Consent for publication

Not applicable.

Competing interests

The authors declare no competing interests.

Publisher's Note

Springer Nature remains neutral with regard to jurisdictional claims in published maps and institutional affiliations.
==== Refs
References

1. Maris JM Recent advances in neuroblastoma N Engl J Med 2010 362 23 2202 2211 10.1056/NEJMra0804577 20558371
Maris JM. Recent advances in neuroblastoma. N Engl J Med. 2010;362(23):2202–11.20558371
2. Tsubota S Kadomatsu K Origin and initiation mechanisms of neuroblastoma Cell Tissue Res 2018 372 2 211 221 10.1007/s00441-018-2796-z 29445860
Tsubota S, Kadomatsu K. Origin and initiation mechanisms of neuroblastoma. Cell Tissue Res. 2018;372(2):211–21.29445860
3. Yeung V Gabriel M Padhye BD Late effects and treatment related morbidity associated with treatment of neuroblastoma patients in a tertiary paediatric centre Cancer Rep (Hoboken) 2023 6 3 e1738 10.1002/cnr2.1738 36268859
Yeung V, Gabriel M, Padhye BD. Late effects and treatment related morbidity associated with treatment of neuroblastoma patients in a tertiary paediatric centre. Cancer Rep (Hoboken). 2023;6(3): e1738.36268859
4. DuBois SG Macy ME Henderson TO High-risk and relapsed neuroblastoma: toward more cures and better outcomes Am Soc Clin Oncol Educ Book 2022 42 1 13 35522915
DuBois SG, Macy ME, Henderson TO. High-risk and relapsed neuroblastoma: toward more cures and better outcomes. Am Soc Clin Oncol Educ Book. 2022;42:1–13.35522915
5. Li Q Long-term survival of neuroblastoma patients receiving surgery, chemotherapy, and radiotherapy: a propensity score matching study J Clin Med 2023 12 3 754 10.3390/jcm12030754 36769402
Li Q, et al. Long-term survival of neuroblastoma patients receiving surgery, chemotherapy, and radiotherapy: a propensity score matching study. J Clin Med. 2023;12(3):754.36769402
6. Wahba A Wolters R Foster JH Neuroblastoma in the era of precision medicine: a clinical review Cancers (Basel) 2023 15 19 4722 10.3390/cancers15194722 37835416
Wahba A, Wolters R, Foster JH. Neuroblastoma in the era of precision medicine: a clinical review. Cancers (Basel). 2023;15(19):4722.37835416
7. Krystal J Foster JH Treatment of high-risk neuroblastoma Children (Basel) 2023 10 8 1302 37628301
Krystal J, Foster JH. Treatment of high-risk neuroblastoma. Children (Basel). 2023;10(8):1302.37628301
8. Shahrajabian MH, Sun Wl. Survey on multi-omics and multi-omics data analysis integration and application. Curr Pharm Anal. 2023;19(4):267–81.
9. Zhu L Identification of potential biomarkers for pan-cancer diagnosis and prognosis through the integration of large-scale transcriptomic data Front Pharmacol 2022 13 870660 10.3389/fphar.2022.870660 35677427
Zhu L, et al. Identification of potential biomarkers for pan-cancer diagnosis and prognosis through the integration of large-scale transcriptomic data. Front Pharmacol. 2022;13: 870660.35677427
10. Sarhadi VK Armengol G Molecular biomarkers in cancer Biomolecules 2022 12 8 1021 10.3390/biom12081021 35892331
Sarhadi VK, Armengol G. Molecular biomarkers in cancer. Biomolecules. 2022;12(8):1021.35892331
11. Shaikh N A review on computational analysis of big data in breast cancer for predicting potential biomarkers Curr Top Med Chem 2022 22 21 1793 1810 10.2174/1568026622666220907121942 36082858
Shaikh N, et al. A review on computational analysis of big data in breast cancer for predicting potential biomarkers. Curr Top Med Chem. 2022;22(21):1793–810.36082858
12. Demir Karaman E Isik Z Multi-omics data analysis identifies prognostic biomarkers across cancers Med Sci (Basel) 2023 11 3 44 37489460
Demir Karaman E, Isik Z. Multi-omics data analysis identifies prognostic biomarkers across cancers. Med Sci (Basel). 2023;11(3):44.37489460
13. Addanki S Applications of circulating tumor cells and circulating tumor DNA in precision oncology for breast cancers Int J Mol Sci 2022 23 14 7843 10.3390/ijms23147843 35887191
Addanki S, et al. Applications of circulating tumor cells and circulating tumor DNA in precision oncology for breast cancers. Int J Mol Sci. 2022;23(14):7843.35887191
14. Langfelder P Horvath S WGCNA: an R package for weighted correlation network analysis BMC Bioinform 2008 9 559 10.1186/1471-2105-9-559
Langfelder P, Horvath S. WGCNA: an R package for weighted correlation network analysis. BMC Bioinform. 2008;9:559.
15. Dennis G Jr DAVID: database for annotation, visualization, and integrated discovery Genome Biol 2003 4 5 P3 10.1186/gb-2003-4-5-p3 12734009
Dennis G Jr, et al. DAVID: database for annotation, visualization, and integrated discovery. Genome Biol. 2003;4(5):P3.12734009
16. Song Z CHDTEPDB: transcriptome expression profile database and interactive analysis platform for congenital heart disease Congenital Heart Dis 2023 18 6 693 701 10.32604/chd.2024.048081
Song Z, et al. CHDTEPDB: transcriptome expression profile database and interactive analysis platform for congenital heart disease. Congenital Heart Dis. 2023;18(6):693–701.
17. Simon N Regularization paths for Cox’s proportional hazards model via coordinate descent J Stat Softw 2011 39 5 1 13 10.18637/jss.v039.i05 27065756
Simon N, et al. Regularization paths for Cox’s proportional hazards model via coordinate descent. J Stat Softw. 2011;39(5):1–13.27065756
18. Robin X pROC: an open-source package for R and S+ to analyze and compare ROC curves BMC Bioinform 2011 12 77 10.1186/1471-2105-12-77
Robin X, et al. pROC: an open-source package for R and S+ to analyze and compare ROC curves. BMC Bioinform. 2011;12:77.
19. Chen B Profiling tumor infiltrating immune cells with CIBERSORT Methods Mol Biol 2018 1711 243 259 10.1007/978-1-4939-7493-1_12 29344893
Chen B, et al. Profiling tumor infiltrating immune cells with CIBERSORT. Methods Mol Biol. 2018;1711:243–59.29344893
20. Yoshihara K Inferring tumour purity and stromal and immune cell admixture from expression data Nat Commun 2013 4 2612 10.1038/ncomms3612 24113773
Yoshihara K, et al. Inferring tumour purity and stromal and immune cell admixture from expression data. Nat Commun. 2013;4:2612.24113773
21. Maeser D Gruener RF Huang RS oncoPredict: an R package for predicting in vivo or cancer patient drug response and biomarkers from cell line screening data Brief Bioinform 2021 22 6 bbab260 10.1093/bib/bbab260 34260682
Maeser D, Gruener RF, Huang RS. oncoPredict: an R package for predicting in vivo or cancer patient drug response and biomarkers from cell line screening data. Brief Bioinform. 2021;22(6):bbab260.34260682
22. Jiang P Signatures of T cell dysfunction and exclusion predict cancer immunotherapy response Nat Med 2018 24 10 1550 1558 10.1038/s41591-018-0136-1 30127393
Jiang P, et al. Signatures of T cell dysfunction and exclusion predict cancer immunotherapy response. Nat Med. 2018;24(10):1550–8.30127393
23. Zhang Q Cell division cycle 20 (CDC20) drives prostate cancer progression via stabilization of beta-catenin in cancer stem-like cells EBioMedicine 2019 42 397 407 10.1016/j.ebiom.2019.03.032 30904606
Zhang Q, et al. Cell division cycle 20 (CDC20) drives prostate cancer progression via stabilization of beta-catenin in cancer stem-like cells. EBioMedicine. 2019;42:397–407.30904606
24. Miles CE Zhu J Mogilner A Mechanical torque promotes bipolarity of the mitotic spindle through multi-centrosomal clustering Bull Math Biol 2022 84 2 29 10.1007/s11538-021-00985-2 35006409
Miles CE, Zhu J, Mogilner A. Mechanical torque promotes bipolarity of the mitotic spindle through multi-centrosomal clustering. Bull Math Biol. 2022;84(2):29.35006409
25. Chinen T NuMA assemblies organize microtubule asters to establish spindle bipolarity in acentrosomal human cells EMBO J 2020 39 2 e102378 10.15252/embj.2019102378 31782546
Chinen T, et al. NuMA assemblies organize microtubule asters to establish spindle bipolarity in acentrosomal human cells. EMBO J. 2020;39(2): e102378.31782546
26. Sutanto R The oncogene cyclin D1 promotes bipolar spindle integrity under compressive force bioRxiv 2023 21 44
Sutanto R, et al. The oncogene cyclin D1 promotes bipolar spindle integrity under compressive force. bioRxiv. 2023;21:44.
27. Boese A Abstract 2334: inhibition of DHFR modulates anti-tumor immunity Cancer Res 2022 82 12_Supplement 2334 2334 10.1158/1538-7445.AM2022-2334
Boese A, et al. Abstract 2334: inhibition of DHFR modulates anti-tumor immunity. Cancer Res. 2022;82(12_Supplement):2334–2334.
28. Gao Y Identification of novel metabolism-associated subtypes for pancreatic cancer to establish an eighteen-gene risk prediction model Front Cell Dev Biol 2021 9 691161 10.3389/fcell.2021.691161 34447748
Gao Y, et al. Identification of novel metabolism-associated subtypes for pancreatic cancer to establish an eighteen-gene risk prediction model. Front Cell Dev Biol. 2021;9: 691161.34447748
29. Wang P Prognostic values of GMPS, PR, CD40, and p21 in ovarian cancer PeerJ 2019 7 e6301 10.7717/peerj.6301 30701134
Wang P, et al. Prognostic values of GMPS, PR, CD40, and p21 in ovarian cancer. PeerJ. 2019;7: e6301.30701134
30. Yao H Lu F Shao Y The E2F family as potential biomarkers and therapeutic targets in colon cancer PeerJ 2020 8 e8562 10.7717/peerj.8562 32117628
Yao H, Lu F, Shao Y. The E2F family as potential biomarkers and therapeutic targets in colon cancer. PeerJ. 2020;8: e8562.32117628
31. Gao Q E2F3 accelerates the stemness of colon cancer cells by activating the STAT3 pathway Front Oncol 2023 13 1203712 10.3389/fonc.2023.1203712 37456248
Gao Q, et al. E2F3 accelerates the stemness of colon cancer cells by activating the STAT3 pathway. Front Oncol. 2023;13:1203712.37456248
