A total of 721 microarray samples, including two datasets named GSE49710 and E-MTAB-8248 (for short, GSE49710 and EMTAB), both detected on Agilent-020382 Human Custom Microarray 44k, were retrieved from NCBI Gene Expression Omnibus and ArrayExpress. The GSE49710 cohort, part of the SEQC project, portrayed Chinese neuroblastoma atlas, while EMTAB displayed German characteristics. Gene expression matrices accompanying clinical information were downloaded directly for the following analyses. F-score, accuracy, sensitivity and specificity were calculated on the whole GSE49710 cohort. Besides, GSE49711, the RNA-seq result of the same samples from GSE49710, was also fetched for lncRNA-related analysis.

To reduce the biases between the two datasets, we normalized the expression levels by equation (1) since the data should be better limited in 0 to 1 in the neural network.

f_i here indicates the expression of each RNA. f i ′ is designated for the transformed level.

After data normalizing, significant gene expression features were selected by chi-square test which is implemented by ‘chi2’ function in the python package sklearn (https://scikit-learn.org/). Genes whose FDR in the chi-square test was less than 0.05 were filtered. Following this principle, only 172 differentially expressed genes were chosen.

Another common feature selection approach, the Principal Component Analysis (PCA), was used to transfer gene expression matrix into principal components. The cox proportional regression was used to filter the components. We then compared the results of the PCA method with the chi-square method.

After feature selection, a classifier was built to classify genes into different subgroups. Genes with similar biological functions were clustered into the same group. The K-means model in sklearn divided the selected genes into two clusters.

Then a supervised classification model based on deep neural networks was built. Both of two feature groups were used as inputs of our classification. The output of this neural network was a patient’s probability which ranged between 0 and 1. 0 would indicate that the patient is likely to be alive and 1 would indicate that the patient would probably be dead.

The structure of this classifier can be seen in Figure 1 . It consisted of two parts, the encoder and the decoder. For the encoder part, we encoded two different groups of features into two 10 dimensional features by two different two-layer networks.

g 1 ′ g 2 ′ are encoded features. w_ij are weights of the networks. f_{g 1} and f_{g 2} are transformed expressions using formula (1) and b_ij are biases of the networks. The function f indicates the activation function which is a nonlinear part of the encoder. Here, we used the ReLU as this nonlinear function:

To combine these two different encoded features together, we applied the attention mechanism to this model (23). With the difference of simply concatenating different features, this attention mechanism can learn the relationship between them.

This step is illustrated in equation 5. G refers to combined features. sig is a nonlinear function.

After seizing combined features, we input them into the decoder part. The decoder part is also a two-layer neural network.

y' is the output of this classifier.

To train this network, we defined the loss function as equation 8.

α and γ are parameters. In this system we set α to 0.2 and γ to 2. y is the true label of each patient. The patients of each group were uneven, so we used Equation 8, which is called the focal loss and was designed to solve this problem instead of cross entropy loss function (24). To avoid over-fitting, we drop 20% neurons in each layer by the dropout method. GSE49710 was chosen to train the network and EMTAB was to validate our model. For GSE49710, 70% of samples were used to train and the rest were to test.

All of the algorithms mentioned in this subsection were realized by tensorflow 2.2 (https://tensorflow.google.cn/). To optimize this neural network, we applied Adam optimizer and set the learning rate to 0.01 (25). The weights of networks were initialized by glorot uniform distribution (26).

To further evaluate our model, we calculated the accuracy, sensitivity, specificity as well as F-score of our model in two cohorts (27).

ROC curves and AUCs were estimated by pROC package in R to assess the performance of the classifier (28). The time-dependent ROC (tROC) curve and its AUC were estimated by survivalROC R package to introduce survival time into our classifier (29). The curves were plotted by ggplot2 R package (30). Besides, traditional cox regression models, devised by Zhong et al. (31) and De Preter et al. (32), were compared with ours.

In order to correspond patients to those two groups of gene features, we performed K-means clustering on patients. ConsensusClusterPlus was used to determine the best k with parameters (cluster algorithmml: km, distance: Euclidean, replicate time: 1000) (33). The CDF plot and the consensus matrix instructed us to cluster patients into 2 groups. A heatmap showing expression levels of features across samples were made by complexHeatmap R package (34). Survival rates between these two subgroups were measured by the log-rank test in KM curves. R packages survival and survminer were used to fit KM equation and plot the curves (35, 36).

In order to test whether our model was independent of other clinical factors and beneficial for clinicians to identify patients’ conditions, we first applied univariable cox regression on the age, MYCN status, gender, tumor stage, INSS-Risk and our probability score. A multivariable cox regression determined whether a covariate involved was decisive. Forest plots were plotted by forestplot package (https://CRAN.R-project.org/package=forestplot). Decision curve analysis was done by ggDCA (https://cran.r-project.org/web/packages/ggDCA/index.html). The construction and plot of the nomogram which can help clinicians to predict survival were done by rms (https://CRAN.R-project.org/package=rms) and regplot (https://CRAN.R-project.org/package=regplot) R package. Finally, an alluvial diagram was used to visualize the characteristics and disease progressions of each patient. This was achieved by the ggalluvial R package (37).

All biological analyses were done on GSE49710. We executed GO enrichment analysis on the two groups of features by clusterprofiler R package respectively (38). Gene Set Enrichment Analysis (GSEA) was done by GSEA software (Broad Institute, Inc., version 4.0.3) with gene set ‘c5.all.v7.1.symbols.gmt’ and default parameters. String (https://string-db.org/) was used to identify protein-protein interactions (PPI) between the 172 features and Cytoscape software was used to visualize the interaction networks. CytoHubba, a module inside Cytoscape, was carried out to identify the hub genes with 12 algorithms (39, 40). For lncRNAs searching, we extracted 250 lncRNA expressions by sorting lncRNA names in gencode.v34.long_noncoding_RNAs.gtf, a collection of known lncRNAs downloaded from GENCODE (https://www.gencodegenes.org/). Only lncRNAs that owned a standard error > 0.2 could be enrolled in the correlation test between the 18 mRNAs. The cut-off values were: p < 0.05 and |coefficient|>0.5. TarBase v.8 and LncBase Predicted v.2 were used to query for mRNA-miRNA and lncRNA-miRNA pairs respectively (41, 42). TarBase v.8 gave 78 CCNB1-binding miRNAs with filters (Species: Homo Sapiens, Method Type: High-throughput, Regulation Type: DOWN, Validation Type: Direct). LncBase Predicted v.2 provided predicted lncRNA-miRNA pairs with cut-off 0.7. Finally, a competing endogenous RNA (ceRNA) network was constructed using Cytoscape.

Inferred abundances of immune cells and normal tissue cells were calculated by single-sample gene set enrichment analysis (ssGSEA) using GSVA R package (43). Gene sets, also known as the markers of each cell, were collected by Charoentong et al. (44). Univariable cox proportional regression tests were exerted on all cells to reveal prognostic immune cells.

For categorical and continuous data with normal distribution, we applied chi-square tests and student t tests to distinguish the differences between groups. When continuous data was not normal distributed, Wilcoxon sum rank tests and ANOVA were utilized. The Pearson correlation test was used to find linear connections between two groups of observations. A p-value<0.05 was considered statistically significant except for emphasis. To account for multiple-testing, the p-values were adjusted using the Benjamini-Hochberg FDR correction. All statistical analyses were two-tailed and done by Python (Python Software Foundation, version 3.8.2) and R (R Foundation, version 3.7.0).

The overall workflow is shown in Figure 1 . After implementing the chi-square test into each feature and the survival in GSE49710, 172 features were selected ( Figure 2A ). The consensus clustering method determined the best k value as 2 to partition the features based on the expression matrix ( Supplementary Figure 1 ). After that, K-means method was used to cluster the features and patients into 2 groups. Fifty genes were the markers of subgroup 1 of 336 patients (for short, S1), whereas the other 122 genes were the markers of subgroup 2 of 162 patients(S2). Basic characteristics were distributed diversely between the two subgroups except for gender ( Supplementary Table 1 ). It is noteworthy that no MYCN amplification was detected in S1 and 92 were detected in S2 in the GSE49710 cohort while only 1 such case was detected in S1 and 45 in S2 in the EMTAB cohort, suggesting that these 172 features and the corresponding subgroups were MYCN-relevant (chi-square test p < 0.001, Supplementary Table 1 ). Also, these subgroups exhibited significant differences in overall survival and event-free survival (both log-rank test p < 0.001, Figures 2B, C and Supplementary Figure 2 ).

To understand the potential biological functions of the genes in each group, we performed GO-enrichment analyses. Notably, many features in group 1 (F1) are related to the JAK-STAT signaling while features in group 2 (F2) aggregated in the cell migration, bone morphogenesis and ubiquitin-protein transferase activities ( Figures 2D, E ). The JAK-STAT pathway promotes tumor cell proliferation, invasion and immunosuppression through a membrane-nucleus cascade (45). A Previous study has shown that the JAK1/2 inhibitor, AZD1480, could abate neuroblastoma tumor cells growth and extend survivals, suggesting that S1 patients not only maintained better survivals with neuroblastoma, but also might potentially respond to drugs such as AZD1480 to recover (46). Next, the GSEA analysis revealed that S1 showed a higher level of metabolism compared to S2 and S2 developed an intensive immune response ( Supplementary Figure 3 ). This might be attributed to the mild symptoms in S1 where patients kept a normal or slightly elevated metabolism. However, accompanying the progression of tumors in S2, the patients started a fierce immune reaction and finally exhausted. These findings suggest that the subgrouping method could help to understand the molecular pathology underlying the differences in prognoses of neuroblastoma patients.

An encoder-decoder model was then trained on the GSE49710 dataset to predict survivals ( Figure 1 ). Since F1 and F2 contributed unequally to the body responses and outcomes, two neural networks were created for them separately in the encoder. In this encoder, a widely used activation function, the ReLU function in the hidden layer; and a binary classification function, the sigmoid function (or say logistic function), in the output layer were employed. The attention mechanism, inspired by human physiology that people would only concentrate on tasks at hand to improve the efficacy of the encoder-decoder framework with rich information, was used to combine the two encoder parts into the decoder (47). The sigmoid function was also used in the final layer, which outputted survival probabilities. If the probability is less than 0.5, we predicted this patient as alive and vice versa ( Supplementary Table 2 ).

To assess the prediction quality of the overall survival status, we calculated the accuracy, sensitivity, specificity and F-score of our model, which reflected the proportion of correct predictions in all samples, true positives in all positives, true negatives in all negatives and the harmonic mean of precision and recall. In the training set, the accuracy (0.918), sensitivity (0.913) and specificity (0.944) were all greater than 0.9, suggesting that it could efficiently forecast whether a patient would be alive or dead using 172 features. Moreover, a descent efficacy was achieved in the test set (accuracy: 0.852, sensitivity: 0.911, specificity: 0.605), however, F1 score was slightly higher (GSE49710: 0.881, EMTAB: 0.886), indicating that the model was suitable for cohorts of various genetic background. The ROC and tROC curves were then generated which further demonstrate the eminence of the neural network ( Figures 3A–C and Table 1 ). The AUCs of the training set and the test set achieved 0.968 and 0.891 respectively ( Figure 3A ), alluding the robustness of our neural network prediction model. Adding survival time into ROC curves, we found that the 5-year-survival AUCs could be boosted to 0.974 in GSE49710 ( Figure 3B ) and 0.896 in EMTAB ( Figure 3C ), which validated that our neural network could classify patients with high precision. Twenty times of 10-fold cross validation showed the stability and robustness of our neural network architecture, suggesting that our attempts of the attention mechanism would be extended into more datasets.

Next, we compared the performance of alternative methods which varied in either feature selection or model construction with our model. To demonstrate that our feature selection was more closely related to prognoses, a broadly used dimension reduction and feature selection method, PCA, was utilized to select features on GSE49710. Kaiser-Harris Criterion suggests that those principal components whose eigenvalue were more than 1 would be retained. In our study, all components had an eigenvalue greater than 1. Variances explained in each component were similar ( Supplementary Figure 4 ). The top 200 principal components were chosen with cumulative variance percent at 89.532% for further analyses. Since the survival data has not been utilized, a univariable cox regression model was implemented to principal components, resulting in 17 components being selected with p < 0.05. The Consensus cluster determined the best k=value as 3 using K-means clustering ( Supplementary Figure 5 ). However, the 3 subgroups were not significantly different in OS ( Figure 4A , log-rank test p = 0.088) but in EFS ( Figure 4B , log-rank test p = 0.029), which indicated that our chi-square-based feature selection could highlight the hub genes in neuroblastoma and partitioned patients into high and low risk groups.

We then compared the performance of our model with several existing models. In order to reduce the biases between the cohorts used in this study and in the literatures, the expression data was normalized before calculating risk scores. We selected two previously published models as well as the ‘gold standard’, MYCN status, and then performed survival prediction in the GSE49710 and EMTAB cohorts ( Supplementary Table 2 ) (31, 32). Our DL-model generated the highest AUC in 5-year tROC curves in both GSE49710 and EMTAB cohorts ( Figures 4C, D ), indicating that our model outperformed the existing model in survival prediction.

To test whether the DL-based prediction model is widely useful for patients of various background conditions, we implemented the univariable cox regressions among the age, MYCN status, gender, diagnostic stage, risk and our output probability. All variables were converted to binaries in this test ( Supplementary Figure 6 ). Only gender failed in this test as p > 0.05, which was discarded in the multivariable regression. In the multivariable cox model, the probability risk was still significant ( Figure 5A ), indicating that our DL-based model had a broad prognostic ability regardless of clinical covariates.

We further used the decision curve analysis (DCA) to evaluate the net benefit of different models (48). We constructed 3 models: only DL-probabilities, only clinical covariates in multivariable cox analysis as well as a combined model. The combined model achieved the highest net benefit no matter how risk threshold was set ( Supplementary Figure 7 ). The data implied that combining DL-model and clinical information could be profitable for clinicians to diagnose and to predict survivals. Therefore, we build up a nomogram which could help clinicians to predict the potential outcomes of the patients beforehand the medical treatments ( Figure 5B and Supplementary Table 3 ). A C-index 0.889 of the nomogram along with the calibration curve predicted 5-year-survival, demonstrated that this scoring system would be handy and practical in the first-line diagnosis ( Figure 5C ).

The alluvial diagram summarized the samples in our study ( Figure 5D ). 69.07% (67/97) of MYCN-amplified patients would be at stage 4 and this tendency was notable (chi-square test: p < 0.001). Only 21 patients diagnosed INSS low risk were classified into S2 and 35 with high risk into S1, showing that our prognostic subgroups were highly clinical-relevant (chi-square test p < 0.001). Subgroups and probabilities were also correlated (chi-square test p < 0.001). In the alluvial diagram, we observed that if patients had MYCN amplified, whether they were old or young, male or female, most of them would be at stage 4, INSS high risk, Subgroup 2. Whereas MYCN was not amplified, an antithetical conclusion would be drawn.

In order to isolate the hub genes of these 172 features, we first retrieved their PPI on the String website ( Figure 6A ). We then input them into Cytoscape software, and used the cytoHubba module to uncover hub genes by using 12 different algorithms. We summed up the top 5 genes in each algorithm and finally selected 18 genes as pivotal molecules in the network.

Next, we aimed to identify lncRNAs that could participate in the regulation of hub genes. GSE49711 is the RNA-seq of the same sample with GSE49710 and was used to uncover potential lncRNA-mRNA pairs. Subsequently, 5 mRNAs (FBXO17, GNG11, CCNB1, KLF2 and CD9) were highly connected with 17 lncRNAs ( Figure 6B and Supplementary Table 4 ). We noticed that both CCNB1 and CD9 could interact with 4 lncRNAs (MYCNOS, TERC, SNHG1, MIR17HG), along with positive coefficients between CCNB1 and lncRNAs and opposite trends on CD9. CCNB1, an oncogene that controls the cell cycle at G2/M, had been found to be overexpressed in hepatocellular carcinoma and pancreatic cancer (49, 50). CD9 is a tetraspanin involved in cell adhesion, metastasis and inflammation in cancer (51, 52). CCNB1 curtailed and CD9 increased survivals in GSE49710 and R2 (https://hgserver1.amc.nl/cgi-bin/r2/main.cgi), which were consistent with previous reports (both p < 0.001, Figures 6C, D ) (49–52).

The ceRNA theory proposed that lncRNAs and mRNAs competed to interact with shared miRNAs, up-regulating downstream RNAs by impairing miRNA activities (53). We created a CCNB1-associated ceRNA network as described in Methods. ( Supplementary Figure 8 , Supplementary Table 5 ). Of note, the mir-302 family (hsa-mir-302-a, -b, -c and -d) which was highly expressed in embryonic stem cells, was associated with CCNB1 and MIR17HG. This indicated that mir-302 might reduce the proliferation of neuroblastoma as it did in other cancers (54, 55).

Since distinct immune response patterns were observed between the two subgroups in GO-enrichment analysis and GSEA ( Supplementary Figure 3 ), we further analyzed the immune microenvironment in neuroblastoma. The ESTIMATE algorithm was used to infer the purity of the microenvironment by scoring immune and stromal cells (56). Two groups did differ in stromal scores but not in immune scores, suggesting that S1 might preserve more normal stromal cells ( Supplementary Figure 9 ). We used the ssGSEA algorithm to convert gene expression data into relative cell proportions ( Supplementary Figure 10A ). The numbers of the T regulatory cells (Tregs), Natural Killer (NK) cells, Monocytes, MDSCs, Eosinophils and central memory CD4 T cells were up-regulated in S2 compared to S1, while the Memory B cells, Macrophages, Gamma Delta T cells and central memory CD8 T cells exhibited opposite trends. Despite a rich amount of innate cytotoxic NK cells in S2, inhibitory immune cells like Tregs and MDSCs might contribute to deficient cytolytic activities. Since GZMA*PRF1 could represent tumor microenvironment cytolytic activities (57), these data indicated that S1 may possess superior cytolytic activities which might eliminate tumor cells conspicuously ( Supplementary Figure 11A ). In addition, only activated CD8+ T cells were connected with survival events in two subgroups at p < 0.05, however, they anticipated contradicting outcomes ( Supplementary Figure 10B , Supplementary Table 6 ). This implied that CD8+ T cells might play dual roles in neuroblastoma patients, i.e., CD8+ T cells functioned as a normal beneficial factor in malignant tumors in S2, however, impeding patients of S1 from recovering.

Then, we examined the intrinsic immune escape mechanism in neuroblastoma. Down-regulations of interferon signals and droppings of two γ-IFN receptors were observed in S2 ( Supplementary Figures 3 , 11B ), whereas IL-2 was increased in S1 which might stimulate T cell differentiation. A loss of HLA-class I/II can aid tumor cells to escape from immune monitoring. HLA-A and HLA-C were lower in S2, making the tumors prone to survive ( Supplementary Figure 11C ) (58, 59). The expression levels of PDCD1, PDL1 and CTLA4, which are critical immune checkpoint genes, were also affected in S2 (60). Overall, these data suggest that disturbance of the immune system may be underneath the poor outcomes of the patients in S2.