Tumor dormancy is closely related to prognosis prediction and tumor immunity in neuroblastoma

Abstract

Background: Neuroblastoma (NB), which is the most frequent and fatal solid tumor in early childhood, lacks an accurate approach to prevent or forecast its recurrence. Dormant NB cells are responsible for metastasis, drug resistance, and suppressive activity in the immune system. However, there is a lack of systematic research on the interaction between dormancy and NB prognosis and its potential associations with tumor immunity.

Methods: We downloaded NB gene expression data and clinical information from the Gene Expression Omnibus and ArrayExpres databases. Based on consensus clustering of the expression of dormancy-associated genes, the NB samples were classified into different groups, and differentially expressed genes (DEGs) were explored in each group. Functional analyses of DEGs were performed, followed by the establishment of a predictive dormancy signature and the assessment of tumor immunity. Finally, sex, age, International Neuroblastoma Staging System (INSS) stage, and MYCN status were identified as independent overall survival-related variables, which were incorporated into the nomogram.

Results: A dormancy-associated gene signature, including CDKN2A, BHLHB3, CDKN2B, MAPK14, CDKN1B, and BMP7, was established. The gene signature showed a strong correlation with NB immune infiltration and capacity to predict NB patient prognosis. A nomogram including MYCN status, INSS stage, age and gene signature risk score was established which further divided NB into high, medium and low-risk groups. This nomogram had certain guiding significance in decision-making for clinical treatment.

Conclusions: Our results suggested that the 6-gene genetic signature for NB based on dormancy could predict NB survival and response to immunotherapy.

Keywords: Dormancy; gene signature; neuroblastoma (NB); prognosis; tumor immunity.

Figure 1

Flowchart showing the research design. KM, Kaplan-Meier; DEG, differentially expressed gene; ROC, receiver operating characteristic.

Figure 2

Neuroblastoma classification based on dormancy-associated genes. (A) Consensus clustering of 498 neuroblastoma samples in the GSE49710 dataset at k=2. (B) Visualization of t-SNE for the 2 clusters. Box plots of the dormancy signature score distribution calculated by ssGSEA (C), PCA (D), and z-score (E) algorithms in the 2 clusters. (F) Kaplan-Meier overall survival curves of the 2 clusters. The log-rank test was used to calculate the P value. Box plots of the DNA replication (G), cell cycle (H), G2/M checkpoint (I), mitotic spindle (J), glycogen synthesis (K), glycolysis (L), and apoptosis (M) score distributions in the 2 clusters. Scores were calculated by the PCA algorithm. *, P<0.05; **, P<0.01; ****, P<0.0001; ns, not significant, two-sided unpaired Wilcoxon test. t-SNE, t-distributed stochastic neighbor embedding; ssGSEA, single-sample gene set enrichment analysis; PCA, principal component analysis.

Figure 3

Immunity analysis in 2 clusters. Distributions of the immune score (A), stromal score (B), and estimate score (C) in the two clusters are shown in box plots. Scores were calculated by the ESTIMATE algorithm. (D) Box plots of immune cell score distribution in the 2 clusters. Scores were calculated by the MCP-Count algorithm. Box plots of the interleukin (E) and interleukin receptor (F) signature score distribution in the 2 clusters. Significantly different immune cells are highlighted in pink. *, P<0.05; **, P<0.01; ***, P<0.001; ns, not significant, two-sided unpaired Wilcoxon test. MCP-Count, microenvironment cell populations-counter.

Figure 4

Construction of predicative dormancy signature. (A) Seven dormancy-associated genes with P<0.01 in the univariate Cox regression analysis of overall survival obtained from the GSE49710 dataset. (B) The forest plot of the dormancy-associated 6-gene signature. (C) Kaplan-Meier overall survival curves of the 6-gene signature. Patients from the GSE49710 dataset were divided into 2 groups based on the median risk score calculated by the equation. The P value was calculated by the log-rank test. (D) Distribution of the risk score, the associated survival data, and the mRNA expression of the 6 genes in the GSE49710 dataset. (E) ROC curves for the 6-gene signature-based predictions of 1-, 3-, and 5-year overall survival. **, P<0.01; ***, P<0.001; ****, P<0.0001. ROC, receiver operating characteristic; AUC, area under the curve.

Figure 5

External validation of predictive signature. (A) t-SNE analysis of the 2 risk-groups derived from the predictive signature in the E-MTAB-8248 dataset. (B) Kaplan-Meier survival curves of the predictive signature. Patients from the E-MTAB-8248 dataset were stratified into 2 groups according to the threshold calculated from the GSE49710 dataset. The log-rank test was utilized to calculate the P value. (C) Distribution of the risk score, the associated survival data, and the mRNA expression of the 6 genes in the E-MTAB-8248 dataset. (D) ROC curves to predict 1-, 3-, and 5-year overall survival based on the 6-gene signature. t-SNE, t-distributed stochastic neighbor embedding; ROC, receiver operating characteristic; AUC, area under the curve.

Figure 6

Clinical relevance and tumor immunity of predictive signature. The risk score distribution among groups with different INSS stage (A), MYCN status (B), and tumor progression status (C) in the GSE49710 dataset. The mRNA expression distributions of NTRK3 (D), LDHA (E), and LDHB (F) in distinct subgroups. Scores were calculated by the ESTIMATE algorithm. Box plots of the distribution of the immune score (G), stromal score (H), and ESTIMATE score (I) in the high- and low-risk groups. Scores were calculated by the ESTIMATE algorithm. (J) Box plots of immune cell score distribution in the 2 clusters. Scores were calculated by the MCP-Count algorithm. *, P<0.05; **, P<0.01; ***, P<0.001 ****P<0.0001; ns, not significant, two-sided unpaired Wilcoxon test. INSS, International Neuroblastoma Staging System; NA, non-amplified; A, amplified; MCP-Count, microenvironment cell populations-counter.

Figure 7

Construction of predictive nomogram. (A) Univariate Cox regression analysis in the GSE49710 dataset. (B) A prognostic nomogram predicting the 1-, 3-, and 5-year overall survival of neuroblastoma. (C) Kaplan-Meier overall survival curves of the nomogram. Base on the nomogram score, patients from the GSE49710 dataset were divided into 3 groups. The P value was calculated by the log-rank test. (D) The distribution of nomogram scores and the corresponding survival information from the GSE49710 dataset. (E) ROC curves for 1-, 3-, and 5-year overall survival predictions for the nomogram model in the GSE49710 dataset. (F) The calibration plot of 1-, 3-, and 5-year survival for internal validation of the nomogram. The Y-axis and X-axis represent the actual overall survival and the predicted overall survival, respectively. ns, not significant. HR, hazard ratio; ROC, receiver operating characteristic; AUC, area under the curve; OS, overall survival.