Multi-omics integration reveals a six-malignant cell maker gene signature for predicting prognosis in high-risk neuroblastoma

Background: Neuroblastoma is the most common extracranial solid tumor of childhood, arising from the sympathetic nervous system. High-risk neuroblastoma (HRNB) remains a major therapeutic challenge with low survival rates despite the intensification of therapy. This study aimed to develop a malignant-cell marker gene signature (MMGS) that might serve as a prognostic indicator in HRNB patients.

Methods: Multi-omics datasets, including mRNA expression (single-cell and bulk), DNA methylation, and clinical information of HRNB patients, were used to identify prognostic malignant cell marker genes. MMGS was established by univariate Cox analysis, LASSO, and stepwise multivariable Cox regression analysis. Kaplan–Meier (KM) curve and time-dependent receiver operating characteristic curve (tROC) were used to evaluate the prognostic value and performance of MMGS, respectively. MMGS further verified its reliability and accuracy in the independent validation set. Finally, the characteristics of functional enrichment, tumor immune features, and inflammatory activity between different MMGS risk groups were also investigated.

Results: We constructed a prognostic model consisting of six malignant cell maker genes (MAPT, C1QTNF4, MEG3, NPW, RAMP1, and CDT1), which stratified patients into ultra-high-risk (UHR) and common-high-risk (CHR) group. Patients in the UHR group had significantly worse overall survival (OS) than those in the CHR group. MMGS was verified as an independent predictor for the OS of HRNB patients. The area under the curve (AUC) values of MMGS at 1-, 3-, and 5-year were 0.78, 0.693, and 0.618, respectively. Notably, functional enrichment, tumor immune features, and inflammatory activity analyses preliminarily indicated that the poor prognosis in the UHR group might result from the dysregulation of the metabolic process and immunosuppressive microenvironment.

Conclusion: This study established a novel six-malignant cell maker gene prognostic model that can be used to predict the prognosis of HRNB patients, which may provide new insight for the treatment and personalized monitoring of HRNB patients.

Introduction

Neuroblastoma (NB) is one of the most common malignant solid tumors of childhood, accounting for 15% of childhood cancer-related deaths (Maris et al., 2007). NB is characterized by extensive heterogeneous clinical phenotype, ranging from spontaneous regression to metastatic disease with poor prognosis (George et al., 2020). The current risk classification system uses clinical and biological variables, such as age, stage, and MYCN oncogene amplification, to stratify NB patients into three risk groups (low-, intermediate-, and high-risk) (Cohn et al., 2009). The survival probability for low-risk patients exceeds 90%, however, it remains below 50% (Ladenstein et al., 2017; Amoroso et al., 2018) in high-risk neuroblastoma (HRNB) patients despite intensive and multi-modal therapy. In addition, patients with an inherently good prognosis but classified into the high-risk group under the current stratification system will undergo toxic treatment, which exposes them to an unnecessary risk of potential long-term side effects (Vermeulen et al., 2009). Therefore, it is imperative to develop more precise biomarkers for HRNB patients to avoid under- or over-treatment and to discriminate who will benefit from new experimental therapy.

The advent of high-throughput sequencing technologies has led to increased efforts to identify molecular prognostic markers in NB from various omics (De Preter et al., 2011; Decock et al., 2012; Garcia et al., 2012; Valentijn et al., 2012; He et al., 2020; Wang et al., 2020), partly with the desire to refine existing risk stratification further. Previous studies have only used dysregulated genes (Chen et al., 2016; Wei et al., 2018) or genomic alterations (Bilke et al., 2008; Stigliani et al., 2012; Depuydt et al., 2018; Fernandez-Blanco et al., 2021) to predict HRNB survival, rarely through multi-omics integration. Additionally, these studies largely used bulk sequencing technologies, which were restricted to mixed cell populations. Therefore, despite successfully identifying many prognostic genes, these studies heavily ignored tumor heterogeneity. NB is a type of tumor intimately related to the early development and differentiation of neuroendocrine (NE) cells (Nunes-Xavier et al., 2021). The intrinsic and extrinsic features of malignant NE cells can be precisely characterized using single-cell RNA sequencing (scRNA-seq) technology, which enables gene profiling and discovery of oncogenic cellular populations and associated marker genes at single-cell resolution (Tang et al., 2009; Patel et al., 2014). Thus, the lack of prognostic stratification for HRNB based on scRNA-seq data and multi-omics integration analysis prompted us to conduct this study.

In this study, we first utilized scRNA-seq data from HRNB patients to identify marker genes of malignant cells. Subsequently, a malignant-cell marker gene signature (MMGS) was constructed to predict the prognosis of HRNB patients based on multi-omics integration analysis. Furthermore, the prognostic value and performance of MMGS were validated in an independent cohort. Finally, the relationships between MMGS and tumor immune/inflammatory features were investigated.

Materials and methods

Data acquisition and preprocessing

ScRNA-seq data of 11 HRNB samples and cell types annotation file were downloaded from GSE137804 (Dong et al., 2020) and were exploited to identify the marker genes of malignant cell. Meta-program gene sets (1 and 6) highly expressed in HRNB and strongly associated with a poor prognosis were downloaded from Dong et al. (2020). Bulk DNA methylation and bulk mRNA expression datasets of 56 HRNB patients were downloaded from GSE73515 and GSE73517, respectively (Henrich et al., 2016). Sequencing Quality Control (SEQC) (Zhang et al., 2015) and Therapeutically Applicable Research to Generate Effective Treatments (TARGET) (Pugh et al., 2013) HRNB cohorts with bulk RNA-seq profiles and clinical information were downloaded from GSE49711 and TARGET datasets,¹ respectively. Gene mutation profiles of HRNB samples were downloaded from the TARGET project.²

The “Seurat” package (Butler et al., 2018) was employed to analyze and visualize scRNA-seq data. The “FindAllMarkers” function was used to identify the marker genes of malignant cells. Adjusted p-value < 0.05 and log2(fold change) > 0.25 were used as the cutoff threshold values to identify marker genes. The “RunUMAP” function was used to embed the cells in a two-dimensional map. For bulk RNA expression data, quantile algorithm from “limma” package (Ritchie et al., 2015) was used to perform normalization. For bulk DNA methylation data, subset quantile normalization was applied according to Touleimat and Tost (2012). We removed the sites missing in more than 30% of the samples and filled the missing values using the k-Nearest Neighbor algorithm. Pearson’s correlation coefficient of each gene corresponding to bulk mRNA expression and bulk DNA methylation profile was calculated. Genes with an absolute correlation coefficient greater than 0.4 and p-value < 0.001 were identified as methylation correlated (METcor) genes. The overlap between malignant cell marker genes and METcor genes was named MK-METcor genes and used for further study.

Construction and validation of malignant cell marker genes prognostic model

A univariate Cox regression analysis was conducted to assess the prognostic value of MK-METcor genes for overall survival (OS) in SEQC HRNB cohort. P-value < 0.05 was used as the cutoff threshold value to identify as a prognostic gene. Next, least absolute shrinkage and selection operator (LASSO) Cox regression (Tibshirani, 1997) analysis was conducted to evaluate the prognostic genes using the “glmnet” package. Finally, a stepwise multivariate Cox regression analysis was used to optimize the prognostic signatures. A linear combination of the risk coefficient and mRNA expression of the genes was used to construct a risk model. Based on the median risk score, HRNB patients were divided into two risk groups: the ultra-high-risk group (UHR group) and the common-high-risk group (CHR group). The Kaplan–Meier (KM) method and the log-rank test were used to calculate the statistical significance of the differences in survival using the “survminer” package. The area under the curve (AUC) of the time-dependent receiver operating characteristic (tROC) curve was calculated to validate the prognostic power of MMGS using the “survivalROC” package (Heagerty and Zheng, 2005). The prognostic ability of MMGS was validated in an independent TARGET HRNB cohort. Clinical model was constructed based on canonical clinical and molecular variables (including age, sex, and MYCN status) by multivariable cox regression.

Differentially expressed gene identification and enrichment analysis

The “limma” package (Ritchie et al., 2015) was used to identify DEGs between UHR and CHR groups. Adjusted p-value < 0.05 and fold change (FC) > 1 were used as the cutoff threshold values to identify differentially expressed genes (DEGs). Compared with the CHR group, DEGs that were overexpressed in the UHR group were considered “overexpressed DEGs” and those that were underexpressed in the UHR group were considered “underexpressed DEGs.”

Gene ontology (GO) and Kyoto Encyclopedia of Genes and Genomes (KEGG) pathway enrichment analyses were conducted using the “enrichGO” and “enrichKEGG” functions in the “clusterProfiler” (Yu et al., 2012) package, respectively. The enriched terms were filtered with the p-value adjusted by the Benjamini-Hochberg method (adjP-value) < 0.05.

Gene Set Enrichment Analysis (GSEA) was performed using the “gseKEGG” and “GSEA” functions in the “clusterProfiler” package (Yu et al., 2012) with gene sets involved in KEGG pathways and hallmarks, respectively. Hallmark gene sets were download from Molecular Signatures Database (MsigDB) (h.all.v6.2.symbols.gmt) (Liberzon et al., 2015). False discovery rate (FDR) q-value < 0.05 and | NES (net enrichment score)| > 1 were set as the significant thresholds.

Tumor immune and inflammatory features analyses

The ESTIMATE algorithm was employed to assess infiltration levels of the immune and stromal cells by the “estimate” package (Yoshihara et al., 2013). The CIBERSORT algorithm (Newman et al., 2019) was used to dissect the proportion of 22 immune cell types infiltration.

Seven metagene clusters (HCK, Interferon, MHC-I, MHC-II, IgG, LCK, and STAT1) were collected from Rody et al. (2009) that have been widely applied to evaluate the tumor inflammatory activity. We conducted a single-sample GSEA (ssGSEA) analysis to obtain an inflammatory activity score by the “GSVA” package (Hanzelmann et al., 2013).

Statistical analysis

All analyses were conducted with R version 4.0.0³ and its appropriate packages. Univariate and multivariate Cox regression analyses were performed to assess the prognostic value of MMGS and various clinical variables. Wilcoxon test or Kruskal–Wallis test was used to analyze differences between different groups. The chi-square test was conducted to investigate the association between MMGS risk groups, MYCN status, and OS status. The “maftools” package (Mayakonda et al., 2018) was used to visualize the genomic mutation landscape. Heatmap plot was visualized using the “pheatmap” package. P-value < 0.05 was set as a significant threshold.

Results

Identification of malignant cell marker genes expression profiles

The overall flowchart of this study was shown in Supplementary Figure 1. First, we obtained scRNA-seq data of 160,910 cells from 11 HRNB patients, including 96,843 malignant NE cells (Dong et al., 2020; Figure 1A). A total of 3,404 malignant cell marker genes were identified (Figure 1B). We observed that malignant cells highly expressed signature genes of chromaffin cells (such as STMN2, TUBB2B, and MEG3), presenting consistent results reported by Dong et al. (2020) that most cancer cells showed a strong chromaffin-cell-like feature. We found that the malignant cell marker genes were enriched in pathways associated with cell cycle and metabolic features, such as oxidative phosphorylation, carbon metabolism, citrate cycle (TCA cycle), regulation of the cellular amino acid metabolic process, and pyruvate metabolism (Figures 1C,D).

FIGURE 1

Figure 1. Identification of malignant cell marker genes by scRNA-seq analysis. (A) UMAP plot colored by various cell types. (B) Heatmap showing the top five marker genes of each cell type. GO biological process (BP) terms (C) and KEGG pathways (D) enrichment analysis of malignant cell marker genes of HRNB. pDC, plasmacytoid dendritic cell.

Construction of malignant-cell marker gene signature prognostic model

Based on the correlation between bulk mRNA expression and bulk DNA methylation, a total of 829 METCor genes were obtained (Supplementary Figure 2A), mainly located on the CpG island (Supplementary Figure 2B). The 168 MK-METcor genes, overlapping between malignant cell marker genes and METcor genes, were used for subsequent analysis (Supplementary Figure 2C).

To construct a prognostic model based on the 168 MK-METcor genes, we first used 176 HRNB patients from the SEQC dataset as the training set to conduct a univariate Cox regression analysis. We identified 35 MK-METcor genes with prognostic value. Subsequently, LASSO Cox regression analysis was used to shrink the variants to 16 genes (Supplementary Figures 3A,B). Finally, stepwise multivariate Cox regression analysis was performed to optimize genes to include only the six most predictive genes. Therefore, MMGS risk score = (0.279 × RAMP1 expression) + (0.278 × CDT1 expression) + (0.024 × NPW expression) + (-0.073 × MAPT expression) + (-0.101 × MEG3 expression) + (-0.347 × C1QTNF4 expression) (Figure 2A). RAMP1, CDT1, NPW, and MAPT showed significant expression differences between CHR and UHR groups, while MEG3 and C1QTNF4 didn’t show expression differences between UHR and CHR groups (Supplementary Figure 3C). The relative expression of these six signature genes in various cell types showed their tumor cell-specific expression (Supplementary Figure 3D). Except for MAPT, the other genes showed a negative correlation between mRNA expression and DNA methylation (Supplementary Figure 3E). By ranking the risk score from low to high, the median risk score (0.901) was used to divide patients into CHR (n = 88) and UHR (n = 88) groups (Figure 2B). The expression profiles of six signature genes in the training set were shown in Figure 2C. KM analysis demonstrated that UHR patients had a significantly shorter OS than CHR patients (p = 0.012, HR: 3.722 [95% CI: 1.800–7.695]) (Figure 2D). A tROC analysis showed that the AUC values for the 1-, 3-, and 5-year OS of MMGS model were 0.780, 0.693, and 0.618, respectively, which were higher than those of the clinical model (except AUC for 5-year OS) (Figure 2E).

FIGURE 2

Figure 2. Establishment of MMGS in the SEQC HRNB cohort. (A) The coefficients of the identified six malignant cell marker genes. (B) The distribution of risk score and survival status. (C) The expression characteristics of the identified six malignant cell marker genes. (D) KM curves of OS between the UHR and CHR groups. (E) ROC curves of MMGS model (solid lines) and clinical model (dot lines) to predict the 1-, 3-, and 5-year OS. ClinModel, clinical model; HR, hazard ratio; CI, confidence interval.

Validation of the prognostic value of malignant-cell marker gene signature

In an independent cohort including 123 HRNB patients from the TARGET project, a risk score was calculated for each patient. Taking the median risk score as the cutoff value, patients were classified into UHR (n = 61) and CHR (n = 62) groups (Figure 3A). KM analysis also showed that the UHR group had an inferior prognosis than the CHR group (Figure 3B, p = 0.033, HR: 3.557, [95% CI: 1.856–6.818]). A tROC analysis showed that the AUC values for 1-, 3-, and 5-year OS of MMGS model were 0.649, 0.669, and 0.681, respectively, which were higher than those of the clinical model (0.552, 0.649, and 0.592, respectively) (Figure 3C).

FIGURE 3

Figure 3. Validation of MMGS in the TARGET HRNB cohort. (A) The distribution of risk score and survival status. (B) KM curves of OS between the UHR and CHR groups. (C) ROC curves of MMGS model (solid lines) and clinical model (dot lines) to predict the 1-, 3-, and 5-year OS. ClinModel, clinical model; HR, hazard ratio; CI, confidence interval.

To further validate the superiority of MMGS, we compared MMGS model to the model constructed without DNA methylation information (not including METcor gene set) (Supplementary Figure 4A). This model consisted of eight signature genes (Supplementary Figure 4B). UHR patients had a significantly shorter OS than CHR patients (p = 0.0011 and 0.042, respectively) in both training and validation sets (Supplementary Figures 4C,D). The AUC values for the 3- and 5-year OS (0.727, 0.675, respectively) of this model were slightly higher than those of MMGS in the training set (Supplementary Figure 4E). However, in the validation set, the AUC values for 1-, 3- and 5-year OS of this model (0.626, 0.632, and 0.562, respectively) were lower than those of the MMGS model (Supplementary Figure 4F).

Additionally, Dong et al. (2020) identified two meta-program gene sets (1 and 6), which were highly expressed in HRNB and strongly associated with a poor prognosis. Then, we compared MMGS model to the model starting with these two meta-program gene sets and with (Supplementary Figure 5) or without DNA methylation information (Supplementary Figure 6). The model constructed with meta-program gene sets and METcor gene set consisted of three signature genes (Supplementary Figures 5A,B). UHR patients had a significantly shorter OS than CHR patients (p = 0.0074) in the training set (Supplementary Figure 5C) and had a nearly significantly shorter OS than CHR patients (p = 0.053) in the validation set (Supplementary Figure 5D). The AUC values for the 3- and 5-year OS (0.72, 0.646, respectively) of this model were slightly higher than those of the MMGS model in the training set (Supplementary Figure 5E). However, in the validation set, the AUC values for 1-, 3- and 5-year OS of this model (0.633, 0.592, and 0.564, respectively) was all lower than those of the MMGS model (Supplementary Figure 5F). The model construction using metaprograms genes without DNA methylation information consisted of eight signature genes (Supplementary Figures 6A,B). UHR patients had a significantly shorter OS than CHR patients (p = 0.00053) in the training set (Supplementary Figure 6C). However, there was no significant difference in OS between UHR and CHR patients (p = 0.81) in the validation set (Supplementary Figure 6D). The AUC values for the 3- and 5-year OS (0.7, 0.706, respectively) of this model were slightly higher than those of MMGS model in the training set (Supplementary Figure 6E). However, in the validation set, the AUC values for 1-, 3- and 5-year OS of this model (0.546, 0.572, and 0.548, respectively) were all lower than those of MMGS model (Supplementary Figure 6F). These results indicated that construction of the model starting with meta-program gene sets had slight better performance in predicting 3- and 5-year OS in the training set. However, it’s performance in the validation set was lower than the MMGS model.

To further investigate whether MMGS model can be used to independently predict the survival of HRNB patients, we conducted univariate and multivariate Cox regression analyses using the risk score, clinical variables, and molecular features. As expected, multivariate Cox regression analysis demonstrated that MMGS model was an independent prognostic factor in the SEQC (Table 1, p = 0.039, HR: 1.585, [95% CI: 1.023–2.457]) and TARGET (Table 2, p = 0.008, HR: 1.910, [95% CI: 1.188–3.070]) cohorts.

TABLE 1

Table 1. Univariate and multivariate analyses of MMGS model in the SEQC HRNB cohort.

TABLE 2

Table 2. Univariate and multivariate analyses of MMGS model in the TARGET HRNB cohort.

The correlation of malignant-cell marker gene signature risk score with clinical variables

To examine the correlation between MMGS risk score and clinical features, we investigated MMGS risk score distribution in patients from the SEQC cohort. We observed that the risk score was lower in the non-MYCN-amplified and alive subgroups whereas higher in MYCN-amplified and dead subgroups (Figure 4A and Supplementary Figure 7A). The sankey diagram was employed to visualize the association between MMGS risk group, MYCN status, and OS status (Figure 4B). We found that patients in the HRNB group had a higher proportion of MYCN amplification with the dead outcome (chi-square test, p = 0.003). Moreover, we found that significantly poorer survival probability was associated with high risk scores in HRNB patients with different clinical features, including age (=5 years), male, and stage 4 (Figures 4C–E and Supplementary Figure 7B). In addition, we also compared the difference in event-free survival (EFS) between UHR and CHR groups and observed almost the same results compared with results of overall survival (Supplementary Figures 7C,D).

FIGURE 4

Figure 4. The relationship between MMGS risk score and clinical variables in the SEQC cohort. (A) The relationship between MMGS risk score and MYCN status, OS status. (B) Sankey diagram showing the association between MMGS risk score, MYCN status, and OS status. KM curves of OS between UHR and CHR patients with age less than 5 years old (C), male (D), and stage 4 (E), respectively. Amp, amplified.

Functional enrichment analysis of malignant-cell marker gene signature risk groups

To explore transcriptomic differences between MMGS risk groups, we further investigated pathways and hallmarks enriched in UHR and CHR groups from the SEQC cohort. Differential expression analysis showed UHR group highly expressed metabolic-related genes, such as HTR1E and KCNH5, while CHR group highly expressed inflammatory-related genes, such as CNR2, TNFSF11, CCL21, and CCL19 (Supplementary Figure 8A). We found that the UHR group significantly enriched in the metabolic-related pathways such as oxidative phosphorylation, and MYC target hallmarks (Figures 5A,C), whereas the CHR group enriched in immune- and inflammatory-related pathways like IL6_JAK_STAT3 signaling, inflammatory response, and interferon alpha/gamma response (Figures 5B,D). Consistent results were also observed in the TARGET cohort (Supplementary Figures 8B–E).

FIGURE 5

Figure 5. Representative pathways enriched in UHR and CHR groups from the SEQC cohort. The top five significantly enriched KEGG pathways in the UHR (A) and CHR (B) HRNB patients. The top five significantly enriched hallmarks in the UHR (C) and CHR (D) HRNB patients.

Malignant-cell marker gene signature risk score was associated with tumor immune and inflammatory features

To assess the potential efficacy of immunotherapies in MMGS risk groups, we explored the relationship between MMGS and tumor immune features in the HRNB microenvironment in the SEQC cohort. Based on the ESTIMATE algorithm, we observed that UHR group had significantly lower immune, stromal, and ESTIMATE scores than CHR group (Figure 6A). Next, using the CIBERSORT algorithm, we found that UHR patients had a higher fraction of T cells follicular helper, NK cells activated, and monocytes but had a lower fraction of B cells naïve, T cells CD4 naïve, and macrophage M1 (Figure 6B). Previous studies have reported that B cells promoted immunotherapy response (Helmink et al., 2020) and macrophages functioned as regulators of tumor immunity and immunotherapy (DeNardo and Ruffell, 2019). Characterization of immune cell infiltration levels might provide implications for the development of immunotherapy targets. Furthermore, we compared the expression levels of immune checkpoints and immune-activity-related genes between UHR and CHR groups. The results demonstrated that most of these genes had relatively lower expression in the UHR group (Figure 6C), suggesting that such tumors might be less responsive to immunotherapies.

FIGURE 6

Figure 6. The association between MMGS and the immune cell infiltration in the SEQC cohort. (A) Differences among stromal score, immune score, and ESTIMATE score between UHR and CHR groups. (B) The comparison of infiltration levels of 22 immune cell types between UHR and CHR groups. (C) Immune checkpoints expression between UHR and CHR groups. *p < 0.05; **p < 0.01; ***p < 0.001; ****p < 0.0001; ns, no significance.

To explore the relationship between MMGS risk scores and tumor inflammatory activities, we investigated the associations between MMGS risk scores and seven metagene clusters. These metagene clusters (HCK, Interferon, MHC-I, MHC-II, IgG, LCK, and STAT1) have previously been reported to represent various tumor inflammatory activities (Rody et al., 2009). We found that most of the metagenes were significantly lower expressed in the UHR group (Figure 7A). Next, we found that the UHR group had a lower inflammatory activity score than the CHR group based on ssGSEA analysis (Figure 7B). Moreover, we observed that MMGS risk score showed negative correlations with all clusters (Figure 7C). Similar results were also observed in the TARGET cohort (Supplementary Figures 9, 10).

FIGURE 7

Figure 7. Relationship between MMGS and inflammatory metagenes in the SEQC cohort. (A) Heatmap showing the expression characteristics of seven clusters of metagenes. (B) The ssGSEA score of inflammatory activity in MMGS risk groups. (C) The correlation between MMGS risk scores and metagenes. *p < 0.05; **p < 0.01; ***p < 0.001.

These significant differences in immune and inflammatory features between the UHR and CHR groups inspired us to explore whether these features can improve the performance of MMGS model. Thus, we combined these features into MMGS to construct seven kinds of models. In the training set, AUC values for 1- and 5-year OS of these models were lower than those of MMGS model; AUC values for 3-year OS were close to that of MMGS model (Supplementary Figure 11A). In the validation set, AUC values for 1- and 5-year OS were close to those of MMGS model; AUC values for 3-year OS were lower than those of MMGS model (Supplementary Figure 11B). These results showed that these immune and inflammatory features didn’t effectively seem to improve the performance of MMGS model.

Mutation landscape of malignant-cell marker gene signature risk groups

We further examined the mutation profiles of MMGS risk groups. ALK stood out that it’s the highest mutation frequency (Supplementary Figure 12A). It was proven that ALK gene mutation was a susceptibility factor and a crucial molecular driver of NB (Mosse et al., 2008). In addition, we found that patients with ALK mutation nearly correlated with poor survival compared to the wild-type (WT) patients in the UHR group (Supplementary Figure 12B, p = 0.06).

Discussion

The majority of HRNB is resistant or refractory to the current intensive multimodality therapy, and the survival rates of these patients remain poor (Irwin et al., 2021). ScRNA-seq technology combined with multi-omics data could help to dissect tumor cell heterogeneity and identify potential prognostic biomarkers effectively. Therefore, we first screened HRNB scRNA-seq data to identify 3,404 malignant cell marker genes and found 35 prognostic MK-METcor genes based on multi-omics integration analysis. After that, we developed a six-malignant cell marker cell signature, termed MMGS, which could serve as a novel and independent indicator for the prognosis of HRNB patients. We found that the CHR group had higher levels of tumor-infiltrating immune cells and inflammatory activity. In addition, we discovered that the genes related to immune checkpoints had higher expression in the CHR group than in the UHR group, indicating that immunotherapy might be more appropriate for CHR patients.

Malignant cell marker genes included in our model are: MAPT, C1QTNF4, MEG3, NPW, RAMP1, and CDT1. MAPT is a gene encoding tau protein, which is closely related to maintaining the function of microtubules. It was reported that MAPT might affect tumorigenesis through dysregulation of cell cycle progression, cell mobility, or organelle organization (Papin and Paganetti, 2020). Zaman et al. (2018) reported that MAPT expression was a biomarker for an increased survival rate in pediatric NB. C1QTNF4 belongs to the C1q and TNF-related family. Previous studies suggested its roles in cancer-related inflammation (Li et al., 2011; Luo et al., 2016) and metabolism (Li et al., 2020; Sarver et al., 2020). MEG3 is the first lncRNA discovered to have tumor suppressor function (Beygo et al., 2015). MEG3 expression is controlled by epigenetics, and it is aberrantly CpG-methylated in various tumors (Modali et al., 2015). MEG3 was reported to be significantly downregulated in NB and negatively associated with the international neuroblastoma staging system (INSS) stage (Ye et al., 2020). Moreover, Novak et al. (2020) reported that MEG3 and EZH2 regulated each other through a negative feedback loop and jointly promoted NB progression, suggesting MEG3 may be a potential treatment target for NB. NPW is a gene encoding neuropeptide protein, which can enhance cortisol secretion from adrenal cells. RAMP1 is a chaperone to the amylin and calcitonin-gene-related peptide (CGRP) receptors. It has been reported that RAMP1 is expressed in various tissues, including the adrenal gland (Nagae et al., 2000; Cottrell et al., 2005). CDT1 plays a pivotal role in cell replication and cell cycle regulation (Yang et al., 2019). The destructive role of CDT1 has been demonstrated in tumorigenesis, progression, and chemoresistance in certain tumor types (Karakaidos et al., 2004; Bravou et al., 2005). These reports suggested that genes identified in our MMGS model might play pivotal roles in malignant cells’ biological behaviors and were potential targets for elucidating molecular mechanisms of HRNB.

In this study, MMGS model was proved to be an independent prognostic tool for HRNB tumors in both training and validation datasets. The excellent prognostic value and performance of MMGS motivated us to explore the potential molecular mechanism. We first conducted GSEA analysis, including KEGG pathway and hallmarks, in MMGS risk groups and observed that the UHR group was significantly enriched in metabolic-related pathways, whereas the CHR group was enriched in immune- and inflammatory-related pathways. Thus, the inferior prognosis of UHR group may be partly due to the dysregulation of the metabolic process, which plays a crucial role in tumor progression (Faubert et al., 2020). Additionally, tumor-infiltrating immune cells have been reported to affect tumor development and prognosis (Petitprez et al., 2020; Paijens et al., 2021). Then, we applied ESTIMATE and CIBERSORT algorithms to compare the levels of immune cell infiltration between UHR and CHR groups. The results demonstrated that the UHR group had a lower level of immune cell infiltration, suggesting that this group has a characteristic of “cold tumors” with reduced antitumor activity (Bonaventura et al., 2019). Low levels of tumor-infiltrating immune cells can facilitate tumor cells to evade immune surveillance and promote malignant progression (Song et al., 2022), which may partially be responsible for the inferior prognosis in UHR patients. Furthermore, we found that most of immune checkpoint genes had relatively lower expression in the UHR group. Besides, we observed that most of the inflammatory-related metagenes also had lower expression in the UHR group. Additionally, MMGS risk score showed negative correlations with all metagenes clusters. These results indicated that UHR patients might be in a more immunosuppressive microenvironment, which may partly account for the decreased OS of UHR tumors. Based on the above findings, we indicated that the underlying mechanism for the powerful prognostic performance of MMGS may be due to the dysregulation of the metabolic process and immunosuppressive microenvironment.

Nevertheless, our study still has some limitations. Our prognostic model was constructed and validated based on retrospective cohorts. A large local cohort is need to further validate our model. Additionally, our study did not analyze the interaction of malignant cell marker genes with microenvironmental cell marker genes. Taken together, although there were still some limitations, our study provided a novel and independent indicator for the prognosis of HRNB patients and illuminated potential molecular mechanisms for HRNB.

Conclusion

In conclusion, a six-malignant-cell marker gene signature was established and validated to have an independent performance in predicting the prognosis of HRNB patients. It may serve as a prognostic signature to better predict the prognosis of HRNB patients and help clinicians to develop personalized treatment for HRNB patients.

Data availability statement

The datasets presented in this study can be found in online repositories. The names of the repositories and accession numbers can be found in the article/Supplementary material.

Author contributions

LZ and JL guided and supervised this work. ZY designed the work and participated in data collecting, data processing, program implementation, and manuscript writing. QL contributed to immune cell infiltration analysis. ZC, JW, and HZ contributed to the article polish. All authors provided critical advice for the final manuscript.

Funding

This work was supported by the Natural Science Foundation of China (NSFC 82070167 and NSFC 81870126).

Acknowledgments

We acknowledge the contributions from the GEO (GSE73515 and GSE137804), SEQC, and TARGET databases.

Conflict of interest

The authors declare that the research was conducted in the absence of any commercial or financial relationships that could be construed as a potential conflict of interest.

Publisher’s note

All claims expressed in this article are solely those of the authors and do not necessarily represent those of their affiliated organizations, or those of the publisher, the editors and the reviewers. Any product that may be evaluated in this article, or claim that may be made by its manufacturer, is not guaranteed or endorsed by the publisher.

Supplementary material

The Supplementary Material for this article can be found online at: https://www.frontiersin.org/articles/10.3389/fninf.2022.1034793/full#supplementary-material

Footnotes

1. ^ https://ocg.cancer.gov/programs/target/data-matrix
2. ^ https://portal.gdc.cancer.gov/
3. ^ http://www.R-project.org

References

Amoroso, L., Erminio, G., Makin, G., Pearson, A. D. J., Brock, P., Valteau-Couanet, D., et al. (2018). Topotecan-vincristine-doxorubicin in stage 4 high-risk neuroblastoma patients failing to achieve a complete metastatic response to rapid cojec: A siopen study. Cancer Res. Treat. 50, 148–155. doi: 10.4143/crt.2016.511

PubMed Abstract | CrossRef Full Text | Google Scholar

Beygo, J., Elbracht, M., de Groot, K., Begemann, M., Kanber, D., Platzer, K., et al. (2015). Novel deletions affecting the Meg3-Dmr provide further evidence for a hierarchical regulation of imprinting in 14q32. Eur. J. Hum. Genet. 23, 180–188. doi: 10.1038/ejhg.2014.72

PubMed Abstract | CrossRef Full Text | Google Scholar

Bilke, S., Chen, Q. R., Wei, J. S., and Khan, J. (2008). Whole chromosome alterations predict survival in high-risk neuroblastoma without Mycn amplification. Clin. Cancer Res. 14, 5540–5547. doi: 10.1158/1078-0432.CCR-07-4461

PubMed Abstract | CrossRef Full Text | Google Scholar

Bonaventura, P., Shekarian, T., Alcazer, V., Valladeau-Guilemond, J., Valsesia-Wittmann, S., Amigorena, S., et al. (2019). Cold tumors: A therapeutic challenge for immunotherapy. Front. Immunol. 10:168. doi: 10.3389/fimmu.2019.00168

PubMed Abstract | CrossRef Full Text | Google Scholar

Bravou, V., Nishitani, H., Song, S. Y., Taraviras, S., and Varakis, J. (2005). Expression of the licensing factors, Cdt1 and geminin, in human colon cancer. Int. J. Oncol. 27, 1511–1518.

Google Scholar

Butler, A., Hoffman, P., Smibert, P., Papalexi, E., and Satija, R. (2018). Integrating single-cell transcriptomic data across different conditions, technologies, and species. Nat. Biotechnol. 36, 411–420. doi: 10.1038/nbt.4096

PubMed Abstract | CrossRef Full Text | Google Scholar

Chen, G., Yang, J., Chen, J., Song, Y., Cao, R., Shi, T., et al. (2016). Identifying and annotating human bifunctional rnas reveals their versatile functions. Sci. China Life Sci. 59, 981–992. doi: 10.1007/s11427-016-0054-1

PubMed Abstract | CrossRef Full Text | Google Scholar

Cohn, S. L., Pearson, A. D., London, W. B., Monclair, T., Ambros, P. F., Brodeur, G. M., et al. (2009). The international neuroblastoma risk group (Inrg) classification system: An inrg task force report. J. Clin. Oncol. 27, 289–297. doi: 10.1200/JCO.2008.16.6785

PubMed Abstract | CrossRef Full Text | Google Scholar

Cottrell, G. S., Roosterman, D., Marvizon, J. C., Song, B., Wick, E., Pikios, S., et al. (2005). Localization of calcitonin receptor-like receptor and receptor activity modifying protein 1 in enteric neurons. Dorsal root ganglia, and the spinal cord of the rat. J. Comp. Neurol. 490, 239–255. doi: 10.1002/cne.20669

PubMed Abstract | CrossRef Full Text | Google Scholar

De Preter, K., Mestdagh, P., Vermeulen, J., Zeka, F., Naranjo, A., Bray, I., et al. (2011). Mirna expression profiling enables risk stratification in archived and fresh neuroblastoma tumor samples. Clin. Cancer Res. 17, 7684–7692. doi: 10.1158/1078-0432.CCR-11-0610

PubMed Abstract | CrossRef Full Text | Google Scholar

Decock, A., Ongenaert, M., Hoebeeck, J., De Preter, K., Van Peer, G., Van Criekinge, W., et al. (2012). Genome-wide promoter methylation analysis in neuroblastoma identifies prognostic methylation biomarkers. Genome Biol. 13:R95. doi: 10.1186/gb-2012-13-10-r95

PubMed Abstract | CrossRef Full Text | Google Scholar

DeNardo, D. G., and Ruffell, B. (2019). Macrophages as regulators of tumour immunity and immunotherapy. Nat. Rev. Immunol. 19, 369–382. doi: 10.1038/s41577-019-0127-6

PubMed Abstract | CrossRef Full Text | Google Scholar

Depuydt, P., Boeva, V., Hocking, T. D., Cannoodt, R., Ambros, I. M., Ambros, P. F., et al. (2018). Genomic amplifications and distal 6q loss: Novel markers for poor survival in high-risk neuroblastoma patients. J. Natl. Cancer Inst. 110, 1084–1093. doi: 10.1093/jnci/djy022

PubMed Abstract | CrossRef Full Text | Google Scholar

Dong, R., Yang, R., Zhan, Y., Lai, H. D., Ye, C. J., Yao, X. Y., et al. (2020). Single-cell characterization of malignant phenotypes and developmental trajectories of adrenal neuroblastoma. Cancer Cell 38, 716–733e6. doi: 10.1016/j.ccell.2020.08.014

PubMed Abstract | CrossRef Full Text | Google Scholar

Faubert, B., Solmonson, A., and DeBerardinis, R. J. (2020). Metabolic reprogramming and cancer progression. Science 368:eaaw5473. doi: 10.1126/science.aaw5473

PubMed Abstract | CrossRef Full Text | Google Scholar

Fernandez-Blanco, B., Berbegall, A. P., Martin-Vano, S., Castel, V., Navarro, S., and Noguera, R. (2021). Imbalance between genomic gain and loss identifies high-risk neuroblastoma patients with worse outcomes. Neoplasia 23, 12–20. doi: 10.1016/j.neo.2020.11.001

PubMed Abstract | CrossRef Full Text | Google Scholar

Garcia, I., Mayol, G., Rios, J., Domenech, G., Cheung, N. K., Oberthuer, A., et al. (2012). A three-gene expression signature model for risk stratification of patients with neuroblastoma. Clin. Cancer Res. 18, 2012–2023. doi: 10.1158/1078-0432.CCR-11-2483

PubMed Abstract | CrossRef Full Text | Google Scholar

George, S. L., Parmar, V., Lorenzi, F., Marshall, L. V., Jamin, Y., Poon, E., et al. (2020). Novel therapeutic strategies targeting telomere maintenance mechanisms in high-risk neuroblastoma. J. Exp. Clin. Cancer Res. 39:78. doi: 10.1186/s13046-020-01582-2

PubMed Abstract | CrossRef Full Text | Google Scholar

Hanzelmann, S., Castelo, R., and Guinney, J. (2013). Gsva: Gene set variation analysis for microarray and Rna-Seq data. BMC Bioinformatics 14:7. doi: 10.1186/1471-2105-14-7

PubMed Abstract | CrossRef Full Text | Google Scholar

He, X., Qin, C., Zhao, Y., Zou, L., Zhao, H., and Cheng, C. (2020). Gene signatures associated with genomic aberrations predict prognosis in neuroblastoma. Cancer Commun. Lond 40, 105–118. doi: 10.1002/cac2.12016

PubMed Abstract | CrossRef Full Text | Google Scholar

Heagerty, P. J., and Zheng, Y. (2005). Survival model predictive accuracy and roc curves. Biometrics 61, 92–105. doi: 10.1111/j.0006-341X.2005.030814.x

PubMed Abstract | CrossRef Full Text | Google Scholar

Helmink, B. A., Reddy, S. M., Gao, J., Zhang, S., Basar, R., Thakur, R., et al. (2020). B cells and tertiary lymphoid structures promote immunotherapy response. Nature 577, 549–555. doi: 10.1038/s41586-019-1922-8

PubMed Abstract | CrossRef Full Text | Google Scholar

Henrich, K. O., Bender, S., Saadati, M., Dreidax, D., Gartlgruber, M., Shao, C., et al. (2016). Integrative genome-scale analysis identifies epigenetic mechanisms of transcriptional deregulation in unfavorable neuroblastomas. Cancer Res. 76, 5523–5537. doi: 10.1158/0008-5472.CAN-15-2507

PubMed Abstract | CrossRef Full Text | Google Scholar

Irwin, M. S., Naranjo, A., Zhang, F. F., Cohn, S. L., London, W. B., Gastier-Foster, J. M., et al. (2021). Revised neuroblastoma risk classification system: A report from the children’s oncology group. J. Clin. Oncol. 39, 3229–3241. doi: 10.1200/JCO.21.00278

PubMed Abstract | CrossRef Full Text | Google Scholar

Karakaidos, P., Taraviras, S., Vassiliou, L. V., Zacharatos, P., Kastrinakis, N. G., Kougiou, D., et al. (2004). Overexpression of the replication licensing regulators Hcdt1 and Hcdc6 characterizes a subset of non-small-cell lung carcinomas: Synergistic effect with mutant P53 on tumor growth and chromosomal instability–evidence of E2f-1 transcriptional control over Hcdt1. Am. J. Pathol. 165, 1351–1365. doi: 10.1016/S0002-9440(10)63393-7

CrossRef Full Text | Google Scholar

Ladenstein, R., Potschger, U., Pearson, A. D. J., Brock, P., Luksch, R., Castel, V., et al. (2017). Busulfan and melphalan versus carboplatin, etoposide, and melphalan as high-dose chemotherapy for high-risk neuroblastoma (Hr-Nbl1/Siopen): An international, randomised, multi-arm, open-label, phase 3 trial. Lancet Oncol. 18, 500–514. doi: 10.1016/S1470-2045(17)30070-0

CrossRef Full Text | Google Scholar

Li, Q., Wang, L., Tan, W., Peng, Z., Luo, Y., Zhang, Y., et al. (2011). Identification of C1qtnf-related protein 4 as a potential cytokine that stimulates the stat3 and Nf-Kappab pathways and promotes cell survival in human cancer cells. Cancer Lett. 308, 203–214. doi: 10.1016/j.canlet.2011.05.005

PubMed Abstract | CrossRef Full Text | Google Scholar

Li, Y., Ye, L., Jia, G., Chen, H., Yu, L., and Wu, D. (2020). C1q/Tnf-related protein 4 induces signal transducer and activator of transcription 3 pathway and modulates food intake. Neuroscience 429, 1–9. doi: 10.1016/j.neuroscience.2019.12.039

PubMed Abstract | CrossRef Full Text | Google Scholar

Liberzon, A., Birger, C., Thorvaldsdottir, H., Ghandi, M., Mesirov, J. P., and Tamayo, P. (2015). The molecular signatures database (Msigdb) hallmark gene set collection. Cell Syst. 1, 417–425. doi: 10.1016/j.cels.2015.12.004

PubMed Abstract | CrossRef Full Text | Google Scholar

Luo, Y., Wu, X., Ma, Z., Tan, W., Wang, L., Na, D., et al. (2016). Expression of the novel adipokine C1qtnf-related protein 4 (Ctrp4) suppresses colitis and colitis-associated colorectal cancer in mice. Cell Mol. Immunol. 13, 688–699. doi: 10.1038/cmi.2016.16

PubMed Abstract | CrossRef Full Text | Google Scholar

Maris, J. M., Hogarty, M. D., Bagatell, R., and Cohn, S. L. (2007). Neuroblastoma. Lancet 369, 2106–2120. doi: 10.1016/S0140-6736(07)60983-0

CrossRef Full Text | Google Scholar

Mayakonda, A., Lin, D. C., Assenov, Y., Plass, C., and Koeffler, H. P. (2018). Maftools: Efficient and comprehensive analysis of somatic variants in cancer. Genome Res. 28, 1747–1756. doi: 10.1101/gr.239244.118

PubMed Abstract | CrossRef Full Text | Google Scholar

Modali, S. D., Parekh, V. I., Kebebew, E., and Agarwal, S. K. (2015). Epigenetic regulation of the lncrna Meg3 and its target C-Met in pancreatic neuroendocrine tumors. Mol. Endocrinol. 29, 224–237. doi: 10.1210/me.2014-1304

PubMed Abstract | CrossRef Full Text | Google Scholar

Mosse, Y. P., Laudenslager, M., Longo, L., Cole, K. A., Wood, A., Attiyeh, E. F., et al. (2008). Identification of alk as a major familial neuroblastoma predisposition gene. Nature 455, 930–935. doi: 10.1038/nature07261

PubMed Abstract | CrossRef Full Text | Google Scholar

Nagae, T., Mukoyama, M., Sugawara, A., Mori, K., Yahata, K., Kasahara, M., et al. (2000). Rat receptor-activity-modifying proteins (Ramps) for adrenomedullin/cgrp receptor: Cloning and upregulation in obstructive nephropathy. Biochem. Biophys. Res. Commun. 270, 89–93. doi: 10.1006/bbrc.2000.2390

PubMed Abstract | CrossRef Full Text | Google Scholar

Newman, A. M., Steen, C. B., Liu, C. L., Gentles, A. J., Chaudhuri, A. A., Scherer, F., et al. (2019). Determining cell type abundance and expression from bulk tissues with digital cytometry. Nat. Biotechnol. 37, 773–782. doi: 10.1038/s41587-019-0114-2

PubMed Abstract | CrossRef Full Text | Google Scholar

Novak, E. M., Gimenez, T. M., Neves, N. H., Vince, C. C., Krepischi, A. C. V., Lapa, R. M., et al. (2020). Meg3 and Meg8 aberrant methylation in an infant with neuroblastoma. Pediatr. Blood Cancer 67:e28328. doi: 10.1002/pbc.28328

PubMed Abstract | CrossRef Full Text | Google Scholar

Nunes-Xavier, C. E., Zaldumbide, L., Mosteiro, L., Lopez-Almaraz, R. Andoin, N, Aguirre, P., et al. (2021). Garcia de protein tyrosine phosphatases in neuroblastoma: Emerging roles as biomarkers and therapeutic targets. Front. Cell Dev. Biol. 9:811297. doi: 10.3389/fcell.2021.811297

PubMed Abstract | CrossRef Full Text | Google Scholar

Paijens, S. T., Vledder, A., de Bruyn, M., and Nijman, H. W. (2021). Tumor-infiltrating lymphocytes in the immunotherapy era. Cell Mol. Immunol. 18, 842–859. doi: 10.1038/s41423-020-00565-9

PubMed Abstract | CrossRef Full Text | Google Scholar

Papin, S., and Paganetti, P. (2020). Emerging evidences for an implication of the neurodegeneration-associated protein tau in cancer. Brain Sci. 10:862. doi: 10.3390/brainsci10110862

PubMed Abstract | CrossRef Full Text | Google Scholar

Patel, A. P., Tirosh, I., Trombetta, J. J., Shalek, A. K., Gillespie, S. M., Wakimoto, H., et al. (2014). Single-cell Rna-Seq highlights intratumoral heterogeneity in primary glioblastoma. Science 344, 1396–1401. doi: 10.1126/science.1254257

PubMed Abstract | CrossRef Full Text | Google Scholar

Petitprez, F., Meylan, M., de Reynies, A., Sautes-Fridman, C., and Fridman, W. H. (2020). The tumor microenvironment in the response to immune checkpoint blockade therapies. Front. Immunol. 11:784. doi: 10.3389/fimmu.2020.00784

PubMed Abstract | CrossRef Full Text | Google Scholar

Pugh, T. J., Morozova, O., Attiyeh, E. F., Asgharzadeh, S., Wei, J. S., Auclair, D., et al. (2013). The genetic landscape of high-risk neuroblastoma. Nat. Genet. 45, 279–284. doi: 10.1038/ng.2529

PubMed Abstract | CrossRef Full Text | Google Scholar

Ritchie, M. E., Phipson, B., Wu, D., Hu, Y., Law, C. W., Shi, W., et al. (2015). Limma powers differential expression analyses for rna-sequencing and microarray studies. Nucleic Acids Res. 43:e47. doi: 10.1093/nar/gkv007

PubMed Abstract | CrossRef Full Text | Google Scholar

Rody, A., Holtrich, U., Pusztai, L., Liedtke, C., Gaetje, R., Ruckhaeberle, E., et al. (2009). T-cell metagene predicts a favorable prognosis in estrogen receptor-negative and Her2-positive breast cancers. Breast Cancer Res. 11:R15. doi: 10.1186/bcr2234

PubMed Abstract | CrossRef Full Text | Google Scholar

Sarver, D. C., Stewart, A. N., Rodriguez, S., Little, H. C., Aja, S., and Wong, G. W. (2020). Loss of Ctrp4 alters adiposity and food intake behaviors in obese mice. Am. J. Physiol. Endocrinol. Metab. 319, E1084–E1100. doi: 10.1152/ajpendo.00448.2020

PubMed Abstract | CrossRef Full Text | Google Scholar

Song, P., Li, W., Guo, L., Ying, J., Gao, S., and He, J. (2022). Identification and validation of a novel signature based on Nk cell marker genes to predict prognosis and immunotherapy response in lung adenocarcinoma by integrated analysis of single-cell and bulk Rna-sequencing. Front. Immunol. 13:850745. doi: 10.3389/fimmu.2022.850745

PubMed Abstract | CrossRef Full Text | Google Scholar

Stigliani, S., Coco, S., Moretti, S., Oberthuer, A., Fischer, M., Theissen, J., et al. (2012). High genomic instability predicts survival in metastatic high-risk neuroblastoma. Neoplasia 14, 823–832. doi: 10.1593/neo.121114

PubMed Abstract | CrossRef Full Text | Google Scholar

Tang, F., Barbacioru, C., Wang, Y., Nordman, E., Lee, C., Xu, N., et al. (2009). Mrna-seq whole-transcriptome analysis of a single cell. Nat. Methods 6, 377–382. doi: 10.1038/nmeth.1315

PubMed Abstract | CrossRef Full Text | Google Scholar

Tibshirani, R. (1997). The lasso method for variable selection in the cox model. Stat. Med. 16, 385–395. doi: 10.1002/(sici)1097-0258(19970228)16:4<385:aid-sim380<3.0.co;2-3

CrossRef Full Text | Google Scholar

Touleimat, N., and Tost, J. (2012). Complete pipeline for infinium((R)) human methylation 450k beadchip data processing using subset quantile normalization for accurate DNA methylation estimation. Epigenomics 4, 325–341. doi: 10.2217/epi.12.21

PubMed Abstract | CrossRef Full Text | Google Scholar

Valentijn, L. J., Koster, J., Haneveld, F., Aissa, R. A., van Sluis, P., Broekmans, M. E., et al. (2012). Functional mycn signature predicts outcome of neuroblastoma irrespective of mycn amplification. Proc. Natl. Acad. Sci. U.S.A. 109, 19190–19195. doi: 10.1073/pnas.1208215109

PubMed Abstract | CrossRef Full Text | Google Scholar

Vermeulen, J., De Preter, K., Naranjo, A., Vercruysse, L., Van Roy, N., Hellemans, J., et al. (2009). Predicting outcomes for children with neuroblastoma using a multigene-expression signature: A retrospective Siopen/Cog/Gpoh study. Lancet Oncol. 10, 663–671. doi: 10.1016/S1470-2045(09)70154-8

CrossRef Full Text | Google Scholar

Wang, Z., Cheng, H., Xu, H., Yu, X., and Sui, D. A. (2020). Five-gene signature derived from M6a regulators to improve prognosis prediction of neuroblastoma. Cancer Biomark. 28, 275–284. doi: 10.3233/CBM-191196

PubMed Abstract | CrossRef Full Text | Google Scholar

Wei, J. S., Kuznetsov, I. B., Zhang, S., Song, Y. K., Asgharzadeh, S., Sindiri, S., et al. (2018). Clinically relevant cytotoxic immune cell signatures and clonal expansion of T-Cell receptors in high-risk mycn-not-amplified human neuroblastoma. Clin. Cancer Res. 24, 5673–5684. doi: 10.1158/1078-0432.CCR-18-0599

PubMed Abstract | CrossRef Full Text | Google Scholar

Yang, C. C., Kato, H., Shindo, M., and Masai, H. (2019). Cdc7 activates replication checkpoint by phosphorylating the Chk1-binding domain of claspin in human cells. Elife 8:e50796. doi: 10.7554/eLife.50796

PubMed Abstract | CrossRef Full Text | Google Scholar

Ye, M., Lu, H., Tang, W., Jing, T., Chen, S., Wei, M., et al. (2020). Downregulation of Meg3 promotes neuroblastoma development through Foxo1-mediated autophagy and mtor-mediated epithelial-mesenchymal transition. Int. J. Biol. Sci. 16, 3050–3061. doi: 10.7150/ijbs.48126

PubMed Abstract | CrossRef Full Text | Google Scholar

Yoshihara, K., Shahmoradgoli, M., Martinez, E., Vegesna, R., Kim, H., Torres-Garcia, W., et al. (2013). Inferring tumour purity and stromal and immune cell admixture from expression data. Nat. Commun. 4:2612. doi: 10.1038/ncomms3612

PubMed Abstract | CrossRef Full Text | Google Scholar

Yu, G., Wang, L. G., Han, Y., and He, Q. Y. (2012). Clusterprofiler: An R package for comparing biological themes among gene clusters. Omics 16, 284–287. doi: 10.1089/omi.2011.0118

PubMed Abstract | CrossRef Full Text | Google Scholar

Zaman, S., Chobrutskiy, B. I., and Blanck, G. (2018). Mapt (Tau) expression is a biomarker for an increased rate of survival in pediatric neuroblastoma. Cell Cycle 17, 2474–2483. doi: 10.1080/15384101.2018.1542898

PubMed Abstract | CrossRef Full Text | Google Scholar

Zhang, W., Yu, Y., Hertwig, F., Thierry-Mieg, J., Zhang, W., Thierry-Mieg, D., et al. (2015). Comparison of Rna-Seq and microarray-based models for clinical endpoint prediction. Genome Biol. 16:133. doi: 10.1186/s13059-015-0694-1

PubMed Abstract | CrossRef Full Text | Google Scholar