Prognostic Signature of Immune Genes and Immune-Related LncRNAs in Neuroblastoma: A Study Based on GEO and TARGET Datasets

Background The prognostic value of immune-related genes and lncRNAs in neuroblastoma has not been elucidated, especially in subgroups with different outcomes. This study aimed to explore immune-related prognostic signatures. Materials and Methods Immune-related prognostic genes and lncRNAs were identified by univariate Cox regression analysis in the training set. The top 20 C-index genes and 17 immune-related lncRNAs were included in prognostic model construction, and random forest and the Least Absolute Shrinkage and Selection Operator (LASSO) regression algorithms were employed to select features. The risk score model was constructed and assessed using the Kaplan-Meier plot and the receiver operating characteristic curve. Functional enrichment analysis of the immune-related lncRNAs was conducted using the STRING database. Results In GSE49710, five immune genes (CDK4, PIK3R1, THRA, MAP2K2, and ULBP2) were included in the risk score five genes (RS5_G) signature, and eleven immune-related lncRNAs (LINC00260, FAM13A1OS, AGPAT4-IT1, DUBR, MIAT, TSC22D1-AS1, DANCR, MIR137HG, ERC2-IT1, LINC01184, LINC00667) were brought into risk score LncRNAs (RS_Lnc) signature. Patients were divided into high/low-risk score groups by the median. Overall survival and event/progression-free survival time were shortened in patients with high scores, both in training and validation cohorts. The same results were found in subgroups. In grouping ability assessment, the area under the curves (AUCs) in distinguishing different groups ranged from 0.737 to 0.94, better in discriminating MYCN status and high risk in training cohort (higher than 0.9). Multivariate Cox analysis demonstrated that RS5_G and RS_Lnc were the independent risk factors for overall and event/progression-free survival (all p-values <0.001). Correlation analysis showed that RS5_G and RS_Lnc were negatively associated with aDC, CD8+ T cells, but positively correlated with Th2 cells. Functional enrichment analyzes demonstrated that immune-related lncRNAs are mainly enriched in cancer-related pathways and immune-related pathways. Conclusion We identified the immune-related prognostic signature RS5_G and RS_Lnc. The predicting and grouping ability is close to being even better than those reported in other studies, especially in subgroups. This study provided prognostic signatures that may help clinicians to choose optimal treatment strategies and showed a new insight for NB treatment. These results need further biological experiments and clinical validation.


INTRODUCTION
Neuroblastoma (NB), one of the most common malignant solid tumors in children, accounts for nearly 15% of childhood cancerrelated deaths (1,2). At present, treatment strategy selection for NB is mainly based on risk stratification. Over the past two decades, the International Neuroblastoma Risk Group (INRG) classification system has defined a unified approach for pretreatment assessment, which includes age at diagnosis, histology, the grade of tumor differentiation, MYCN status, 11q aberration, and ploidy (3). According to the classification criteria, low-and intermediate-risk patients can receive reduced chemotherapy and attain good outcomes. However, the prognosis for high-risk patients is still dismal, even though they have received intensive chemotherapy, surgery, radiotherapy, and other combined therapies, and the long-term survival of these patients remains less than 50% (4)(5)(6). Besides, the clinical use of these markers may be limited in some patients despite the elaborate risk stratification.
Thanks to the development of microarrays and highthroughput technologies, an increasing number of abnormalities have been discovered in both driver and nondriver genes in NB (7). Unlike adult cancers, the total somatic mutation rate is less than 20%, even in the high risk NB, as reported by a study from TARGET enrolling 240 patients. Somatic mutations are most commonly seen in ALK (9.2%), followed by PTPN11 (2.9%), ATRX (2.5%), and MYCN (1.7%) (8). Typically, ARID1A and AIRD1B alterations have been identified through integrated genomic analysis, and they are correlated with the poor response to early treatment and unfavorable outcomes (2). Therefore, targeted mutation therapy is not suitable for tumors in children with low mutation frequency, such as NB. In this regard, it is still necessary and urgent to identify new molecular targets and to improve the prognostic outcomes in children with NB.
The tumor microenvironment (TME) has been confirmed to play an important role in various cancer types, including NB (4). The TME components, such as immunocytes (9) and stromal cells, participate in numerous cancer progresses, such as tumorigenesis, metastasis (10,11), response to immune checkpoint blocking therapy (12), and prediction of clinical outcomes (13,14). It is shown that the spontaneous regression in stage 4s NB may be associated with cellular immunity (15). Besides, some researchers have proved that immunocytes like tumor-associated macrophages (TAMs) and cancer-associated fibroblast (CAF) are markedly correlated with the unfavorable consequences in NB (16). Moreover, the impact of tumorinfiltrating immunocytes on the prediction of NB prognosis has also been examined (17,18).
Liu Z depicted the immune gene expression profiles in high risk NB and found an ultra-high risk NB (UHR-NB) group that showed the worst prognosis. Besides, the authors identified four up-regulated immune genes (ADAM22, GAL, KLHL13, and TWISTT1) among the UHR-NB patients (19). However, to date, the prognostic value of immune-related genes and lncRNAs in NB has not been elucidated yet. The anti-GD2 monoclonal antibody is the only immune therapeutic agent used for high risk NB patients, which improves the 2-year overall survival (OS) and event-free survival (EFS). However, it is associated with drawbacks such as significant cumulative toxicity and its improvement of long-term survival is unclear (20). It is noteworthy that the low toxicity of anti-tumor agents and long-term survival are quite important for children.
The present study aimed to explore the prognostic immunerelated genes and lncRNAs, and to discuss their feasibility as prognostic and follow-up markers in NB to help clinicians optimize individual treatment protocols. This study focused on the significance of immune-related prognostic signatures in predicting outcomes for NB patients in each group.
If one gene corresponded to multiple probes, the highest expression value was taken as the representative and log2 transformed. For lncRNAs, we confirmed them as "ncRNAs" from the NCBI database (https://www.ncbi.nlm.nih.gov/gene/).

Identification of Survival-Related Immune Genes
The immune-related gene list was downloaded from the ImmPort database (21) (https://immport.niaid.nih.gov/), including 1534 genes, after excluding duplicates and mismatches, a total of 993 genes remained. Univariate Cox proportional model survival analysis (22) was performed in GSE49710. Genes with p <0.01 were considered as survivalrelated immune genes.

Construction of the Immune-Related Prognostic Signatures
Univariate Cox regression analysis was used to select candidate immune-related genes and rank them according to the Concordance index. The Random Forest algorithm (23,24) was employed to identify the importance of the top 20 C-index immune-related genes in vital status. We then carried out multivariate Cox regression and obtained the coefficient of the genes. At last, a prognostic score model was constructed, including five genes (CDK4, PIK3R1, MAP2K2, ULBP2, and THRA).
where bi represented the coefficient of each gene got from multivariate Cox regression analysis of overall survival time, and xi indicated the expression value by log2 transformed of the corresponding gene.
All patients had risk scores and were divided into high-risk and low-risk score groups according to the median value. The Kaplan-Meier curves of overall survival and event/progression/ relapse-free survival were plotted.

Extraction of Immune-Related LncRNAs and Construction of Prognostic lncRNA Signature
LncRNAs correlated with the five genes which Spearman Correlation Coefficient >0.5 as immune-related lncRNAs. We chose the Spearman correlation because of the non-normal distribution of the gene expression values. The Least Absolute Shrinkage and Selection Operator (LASSO) regression analysis (25) was performed to screen the prognostic immune-related lncRNA signature, and minimal lambda was considered optimal. The modeling process is the same as above.

Functional Annotation and Enrichment Analysis
To explore the biological functions of prognostic immune genes and immune-related lncRNAs, KEGG (26) (Kyoto Encyclopedia of Genes and Genomes) in the STRING (27) database (https://string-db.org/) was performed. FDR <0.05 was considered significant.

Cell Fraction Calculation in Tumor Microenvironment
To describe the cell composition of the tumor microenvironment (TME) in neuroblastoma, we employed the xCell algorithm (28) (http://xcell.ucsf.edu/). The xCell scores of cell types represented a fraction of cells in tumor tissues.

Statistical Analysis
In this study, the R software version 3.6.2 and the GraphPad Prism version 8.0 were used. Univariate and multivariate Cox regression analyses were performed by R package "survival", at the same time, Harrell's concordance index (C-index) was calculated to assess the performance of the model. The ggforest plots were generated by R package "survminer". Kaplan-Meier algorithm was used to evaluate the survival time of high/low-risk groups, and a Log-rank test was employed. To test the predictive ability of the prognostic signature, the area under the curves (AUC) of receiver operating characteristic (ROC) was calculated and plot by R package "pROC". The random forest algorithm was carried out by R package "randomForest". The Lasso regression analysis was implemented by R package "glmnet". The Spearman Correlation Coefficient was calculated by R method "cor.test ()", and the correlation heatmaps were plotted by GraphPad Prism 8.0. The Mann-Whitney U test was used to compare the xCell score of cells in different groups. All statistical tests were two-sided and a p-value < 0.05 was considered statistically significant.

Identification of Survival-Related Immune Genes
Univariate Cox regression analysis was performed in the GSE49710 dataset. Altogether 681 immune-related prognostic genes associated with OS time were identified, and a p-value <0.01 indicated statistical significance. Among these identified genes, 597 were positively correlated whereas 84 were negatively correlated with survival. As suggested by functional enrichment analysis, these genes were mainly enriched in tumor-related pathways (like the PI3K-AKT, MAPK, Ras signaling pathways) and immune-related pathways (such as NK cell-mediated cytotoxicity and T cell receptor signaling pathway), as shown in Supplementary Figure 1. The top 20 prognostic immunerelated genes in terms of their C-indexes were enrolled for subsequent model screening ( Table 1).

Feature Selection by Random Forest
The immune-related genes in Table 1 were included in the prognostic model feature importance selection. Meanwhile, the random forest algorithm was employed to select the important immune-related genes. The above-mentioned 20 genes were used as the features to discriminate the vital status, then the mean decrease accuracy (MDA) and mean decrease Gini (MDG) were obtained. After repeating the algorithm 100 times, the average out-of-bag (OOB) value was 16.99%. As shown in Figure 1A

Prognostic Model RS5_G Had Good Performance in Predicting Neuroblastoma Outcomes
Patients in the GSE49710 dataset (training cohort) were divided into high-or low-risk score groups according to the median RS5_G prognostic score. As shown in Figure 2, the OS time and EFS time of the high risk group were greatly inferior to those of low risk group. To assess the predictive performance, stratification survival analysis was performed in clinical subgroups. In unfavorable clinical subgroups (including stage 4, MYCN amplified, high risk, age ≥ 18 months, progression, and class label unfavorable groups), the OS time of patients with high prognostic scores was significantly worse (all p-values <0.0001).
Moreover, the grouping ability of the prognostic signature was also evaluated by the receiver operating characteristic (ROC) curve analysis.  Figure 3A). To verify the associations between RS5_G (including the immune genes) and clinical characteristics, their corresponding Spearman correlation coefficient (SCC) was calculated, as shown in Figure 3B. The results showed that RS5_G, CDK4, and MAP2K2 are positively correlated to unfavorable prognostic factors, including age group, MYCN amplified, INRG high risk, INSS stage, and progression, which were negatively correlated to OS and EFS days, whereas the expression of PIK3R1, THRA, and ULBP2 showed the opposite results.
The following parameters were incorporated in both univariate and multivariate Cox regression analyses, including RS5_G, age group, high risk, MYCN status, and INSS_h1 (0 for stage 1/2/4s, 1 for stage 3/4). The results showed that the abovementioned features were all risk factors. The _G high, age ≥18 months, high risk, MYCN amplified, and INSS stage 3/4 were correlated with unfavorable outcomes. However, the RS5_G signature had the highest C-index ( Table 2). Upon multivariate Cox analysis, RS5_G and high risk were identified  as the independent prognostic factors for OS, whereas RS5_G, high risk, and INSS_h1 were the independent prognostic factors for EFS. However, RS5_G had the most significant p-value ( Table 3). We then compared the above risk factors with RS5_G seriatim, the results revealed that RS5_G is the only independent prognostic risk factor (Supplementary Table 1).

Immune-Related Prognostic lncRNA Signature Performed Well in Predicting Outcomes and Grouping Different Prognostic Groups of Neuroblastoma
The SCCs were calculated between lncRNAs and genes in the RS5_G signature. The lncRNAs with SCC >0.5 and univariate Survival analysis was subsequently implemented for RS_Lnc, and the results showed that patients with high scores had a significantly shorter OS and EFS time than those with low scores. The results of subgroup analysis showed that the OS time of patients in high risk groups for all subgroups was significantly shortened, except for MYCN amplified groups ( Figure 4).
Thereafter, the ROC curve analysis was adopted to assess performance in discriminating different prognostic groups by RS_Lnc. As a result, the AUC values of age group, high risk,  Figure 5B. The results indicated that RS_Lnc and DANCR are positively correlated to inferior prognostic factors, including age >18 months, MYCN amplified, INRG high risk, advanced stage, and progression, which were negatively correlated to OS and EFS days, whereas the expression of the other ten lncRNAs displayed the opposite results.
Upon univariate Cox analysis, RS_Lnc was associated with short OS and EFS time (HR: 2.718, 1.85; p-value: < 2E-16). Multivariate Cox analysis identified that RS_Lnc and high risk were the independent prognostic factors for OS, whereas RS_Lnc, high risk, and INSS_h1 were the independent prognostic factors for EFS. However, RS_Lnc had the most significant p-value (Supplementary Table 2).

Immune-Related lncRNAs Mainly Enriched in Cancer-and Immune-Related Pathways
To explore the biological functions of immune-related lncRNAs in the prognostic signature, genes with lncRNAs SCCs > 0.5 were screened and the online database STRING was utilized to obtain the protein-protein interaction (PPI) network. The KEGG pathways are shown in Figure 6 and Supplementary Figure 4. These lncRNAs were mainly enriched in some molecular signaling pathways (such as PI3K-AKT, Ras, MAPK, and ErbB signaling pathways) and some immune-related pathways (like T cell receptor signaling pathway and natural killer cellmediated cytotoxicity).

Immune Infiltration and Correlation With Prognostic Signature
In this study, the xCell algorithm was employed to calculate the xCell scores of different cell types in the tumor tissues of NB patients from the GSE49710 dataset, which could represent the proportions of various cells in the tumor microenvironment (TME). Typically, the activated dendritic cells (aDC), T helper cells (Th2 cells), and CD8+ T cells, and their correlations with the prognostic signatures RS5_G and RS_Lnc were analyzed, respectively. These three cell types were chosen because of their higher xCell scores and their associations with tumor prognosis. As shown in Figure 7, aDC and CD8+ T cells were negatively correlated with RS5_G and RS_Lnc, whereas Th2 cells were positively correlated with the above two prognostic signatures. There were significant differences in these three cell types between the high-and low-risk score groups.

Performance Validation in the Independent Datasets
To further confirm the robustness of the prognostic signature RS5_G, its performance was validated in another three independent datasets, GSE16476, TARGET-Asgharzadeh -249, and GSE85047. Furthermore, the performance of RS_Lnc was also validated in GSE16476. The OS and progression-free survival (PFS) time were significantly shorter in high-risk score groups than in low-risk groups of the GSE16476 dataset, both for RS5_G and RS_Lnc  Figures 9A, B). Univariate Cox regression analysis identified that RS5_G and RS_Lnc were the risk factors for OS time in GSE16476, which had higher C-indexes than the other risk factors (age group, MYCN amplified, and INSS_h1) (Supplementary Table 3). Upon multivariate Cox analysis, RS5_G, RS_Lnc, and age group were identified as independent risk factors. However, our signature had a more significant p-value (Supplementary Table 4).
In the TARGET-249 and GSE85047 cohorts, some immunerelated lncRNAs were not detected. Therefore, the robustness of the RS5_G signature was validated. In TARGET-249, two patients with no available survival information were excluded, and finally, 247 patients were included in subsequent analysis. The results suggested that patients with high-risk scores had inferior outcomes to those with low-risk scores. In subgroups analysis including INSS stage 4, event, high risk, age ≥18 months, pathological diagnosis as NB, undifferentiated or poorly differentiated, and unfavorable histology, patients with highrisk scores had visibly reduced OS time (Supplementary Figure 5). In the GSE85047 dataset, a total of 268 patients were included. It was found that OS and PFS were significantly shortened in high-risk score groups (Supplementary Figure 6). Subgroup analysis also came to consistent results. RS5_G had the highest C-indexes in univariate Cox regression analysis of OS and PFS (Supplementary Table 5). Moreover, results of multivariate Cox analysis displayed that, RS5_G and INSS_h1 were the independent risk factors; however, RS5_G had a more significant p-value (Supplementary Table 6

DISCUSSION
Neuroblastoma (NB) is a malignant solid tumor in children, which is linked with highly heterogeneous outcomes ranging from spontaneous regression without treatment to rapid progression after intensive multidisciplinary therapy (3,29).  NB accounts for 8-10% of all malignancies in children and is responsible for 15% of pediatric cancer-related deaths (4,5,19,30). Within the past few decades, the MYCN oncogene has been extensively and deeply studied, and more molecules have also been discovered to play important roles in tumor development, progression, and metastasis, including ALK (8), LMO1 (31, 32), LIN28B (33,34), AURKA (35,36), AIRD1A and ARID1B (2), etc. Besides, some gene inhibitors have been attempted to treat NB, like ALK and AURKA. Nevertheless, the long-term survival of NB at stage 4 or high risk has not been significantly improved yet, and there are still many relapsed or refractory patients. Hence, it is necessary to identify the novel molecules as the therapeutic targets or prognostic biomarkers, especially for high risk patients. At present, NB treatment is mainly based on the INRG risk stratification system (3). Researchers have discovered some prognostic markers that may assist clinicians in choosing the optimal treatment regimen, but few patients benefit from these markers, especially those in subgroups (37)(38)(39)(40). It has not yet been elucidated whether the prognosis of NB patients can be predicted by immune-related genes and immune-related lncRNAs and their correlations with clinicopathological characteristics. In this study, it was found that over 1/8 immune-related genes were positively correlated with NB survival, which meant that the expression of these genes was down-regulated in patients with poor prognosis, indicating that patients with poor prognosis might be in an immunosuppression state. Two immune-related prognostic signatures RS5_G and RS_Lnc were identified, and their predicting and grouping abilities were analyzed in another three independent datasets. The associations of these two signatures with clinicopathological features and some immunocytes were also examined.
Five genes, including CDK4, PIK3R1, THRA, MAP2K2, and ULBP2, were included in the immune-related gene signature. The use of CDK4 inhibitor alone or in combination with other drugs in the treatment of NB has entered a phase I clinical trial (41)(42)(43). MAP2K2, also referred to as MEK2, is an important molecule involved in the MAPK signaling pathway. The MEK inhibitor binimetinib combined with the CDK4/6 inhibitor ribociclib achieved a synergistic therapeutic effect in the preclinical model of high risk NB (42). In this study, CDK4 and MAP2K2 were negatively correlated with long-term survival time. PIK3R1, which encodes the p85a (the regulatory subunit of PI3K), and plays an important role in the PI3K-AKT signaling pathway (44). Fulda elaborated on the feasibility of using the PI3K-AKT-mTOR pathway as a therapeutic target in NB (45). In our study, PIK3R1 was positively correlated with OS and EFS, which is consistent with its function. The above results demonstrated the accuracy and effectiveness of our prognostic gene screening method. THRA is a gene encoding the specific high-affinity nuclear receptor that mediates thyroid hormone. Studies have shown that THRA and THRB (another gene that encodes this receptor), may be involved in human cancer (46). For instance, a study from Xiang Z et al. demonstrated that in ERBB2+ gastric cancer (GC), THRA might be one of the sensitive biomarkers for the targeted ERBB2 treatment, and its expression level was negatively correlated with Myc activation (47). In this study, it was discovered that patients with high THRA expression displayed favorable outcomes. NKG2D is an activating immune receptor found in CD8+ T cells and NK cells, and its ligands include ULBPs, MICA, and MICB. Raffaghello et al. found through experiments, that ULBP2 was expressed in half of primary tumors, and concluded that NB evaded immune cell killing by down-regulating and/or releasing the NKG2D ligands (48). In our study, patients with poor prognosis had reduced ULBP2 expression, and high ULBP2 expression was associated with longer OS and EFS time.
The RS5_G signature constructed in this study performed well in predicting prognosis, both in training and validation cohorts, especially in subgroups. Taking the performance of RS5_G in distinguishing the survival differences in INRG high risk group as an example, our signature achieved comparable results with fewer genes compared with previous studies. Liu identified an ultra-high risk NB (UHR-NB) group by employing the spectral co-clustering algorithm based on 283 immunerelated genes and the two-year survival rates were 48.8% and 43.2% for TARGET and GSE49710, respectively; whereas the two-year survival rates were 79% and 75% in HR-NB group, respectively (19). In our study, the high risk patients were divided into high or low risk score groups according to the median. In the high-risk score groups, the two-, three-and five-year survival rates were 58.7%, 45.0% and 33.9% in TARGET; while 58.0%, 43.2% and 37.5% in GSE49710, respectively. In low-risk score groups, the two-, three-, and five-year survival rates were 78.7%, 64.8%, and 50.9% in TARGET, whereas 83.9%, 78.2%, and 64.4% in GSE49710, respectively. The two-year survival rates of our groups were higher than Liu Z for the unfavorable prognostic groups, yet, only five genes were used in our study, which achieved higher feasibility in clinical trials. Notably, the survival time was distinctly different between high-and lowrisk score groups. Besides, the genes included in our signature are different from those previously published by Liu et al, we considered that it is due to our different methods and concerns. By comparing acute lymphoblastic leukemia (ALL), acute myeloid leukemia (AML), Wilm's tumor (WT), and NB, Liu et al. screened out NB-specific immune-related genes and paid attention to their relationship with prognosis, whereas we focused on immune-related genes that are closely related to NB survival. However, the two signatures screened by us and Liu et al. screen out UHR-NB patients well. Moreover, the robustness of the RS5_G signature was also validated in another two independent cohorts, GSE16476 and GSE85047, and the results were satisfying.
The results of multivariate Cox regression analysis indicated that RS5_G, high risk, and INSS_h1 can be independent prognostic factors for OS and EFS in GSE49710, and RS5_G had the most significant p-value. It demonstrated that the RS5_G signature might help clinicians to choose the appropriate therapeutic regimen for different patients. For example, we can choose a stronger plan at the time of initial treatment for patients with high-risk scores, rather than the intensive treatment at the time of relapse. By comparing RS5_G with the characteristics previously included in NB risk stratification through multivariate Cox analysis, we believe that the immune-related prognostic signature screened by gene expression profiles can help clinicians distinguish patients with different outcomes in INRG high risk groups, to give different treatment regimens. Of course, further experimental verification is needed for accurate application in clinical practice, in particular, the cut-off value is needed to confirm by moderate methods.
Long non-coding RNAs (lncRNAs) participate in the pathological process of tumors including NB in a variety of ways (6,49), and lncRNAs may serve as the prognostic markers for various cancers (37). For example, NBAT1 is located on chromosome 6p22, regulates NB progression, and is correlated with the favorable outcome of NB (50). In addition, LncNB1 plays an oncogenic role in NB tumorigenesis (51). Gao L constructed a risk signature comprised of 10 lncRNAs, with the AUC value of as high as 0.941 in predicting patients' survival (37). However, its generalization performance was unsatisfactory and the signature was not further verified in more independent datasets. In this study, an immune-related lncRNA prognostic signature incorporating 11 lncRNAs was constructed and called RS_Lnc. It served as an independent prognostic factor for OS and EFS in GSE49710 and OS in GSE16476 and had the most significant p-value. Gao identified 10 lncRNAs through differentially expressed analysis and incorporate them in a prognostic signature (37). By contrast, the 11 lncRNAs identified in this study were not the same as the results of the study by Gao, which might be ascribed to the different dataset detection and screening methods.
These eleven lncRNAs were associated with good OS and EFS, except for DNACR that was related to adverse outcomes. Few of them are correlated with cancer tumorigenesis and metastasis. For example, DANCR has been proven to participate in several cancer pathological processes, including tumor cell proliferation, invasion, metastasis, and chemo-resistance, and it may serve as a therapeutic target (52,53). FAM13A1OS (FAM13A-AS1) has been discovered to show high expression in low risk NB patients, which is associated with autophagy and long-term survival (54). Furthermore, DUBR has also been reported to be downregulated in high risk NB and differentially expressed between favorable and unfavorable outcomes (55). MIAT promotes the growth and metastasis of GC and colorectal cancer (CRC) cells via different pathways (56,57), and affects the response of breast cancer (BC) cells to estrogen (58). LINC00667 can promote cell proliferation and migration in non-small cell lung cancer (NSCLC) (59), CRC (60), and Wilm's tumor (61). However, in this study, MIAT and LINC00667 were correlated with favorable outcomes, and the functions of MIAT in NB should be further confirmed in biological experiments. MiR137 plays a role as a tumor suppressor in cancer; meanwhile, HSF1 can directly interact with its host gene MIR137HG to inhibit its expression and promote CRC occurrence (62).
In our study, functional enrichment analysis showed that lncRNA-associated genes are mainly enriched in the cancerrelated signaling pathways (like PI3K-AKT, Ras, MAPK signaling pathways) and immune-related pathways (such as the chemokine and T cell receptor signaling pathways, NK cellmediated cytotoxicity pathway). This indicates that in-depth research from the aspects of signaling pathways and immune response can be considered in the treatment of NB.
Immunocytes, including CD8+ T cells, NK cells, Th2 cells, etc., may assist patients in killing cancer cells or help tumor cells to evade immune killing (11,20). In our previous study, we discovered that CD8+ T cells and aDCs are correlated with favorable outcomes, whereas Th2 cells are associated with unfavorable survival in NB (63). Fruci and colleagues have reported that intratumoral DCs and NK cells are associated with increased T-cell infiltration and favorable survival of NB, and identified two immune gene signatures associated with DC and NK cells that can predict the prognosis of NB patients (64). In this work, the xCell algorithm was employed to mimic immune cell proportion, and the correlations of Th2 cells, CD8+ T cells, and aDCs with the prognostic signature were calculated. It was found that the risk scores were remarkably negatively correlated with CD8+ T cells and aDCs, but positively associated with Th2 cells. These results demonstrated that both innate immunity and acquired immunity may play important roles in NB. The correlations of RS5_G and RS_Lnc with clinicopathological features were examined. As a result, these two signatures were positively correlated with known risk factors of NB but negatively associated with OS and EFS time.
Moreover, the grouping ability of our signature was verified in three datasets, GSE49710, GSE16476, and GSE85047. The AUC values in distinguishing different groups ranged from 0.737 to 0.94. Typically, the signature performed well in discriminating MYCN status and vital status compared with other groups. In the future, further research is warranted to examine the relationship between genes/lncRNAs and MYCN amplified and to explore the possibility of these genes as therapeutic targets.
Even though our prognostic signatures performed well in both the training set and validation cohorts, certain limitations should be noted in this study. Firstly, the datasets came from different platforms, and there were differences in testing, which might cause some bias in the results. Secondly, some lncRNAs were not tested in all datasets, and the performance of the RS_Lnc signature was not verified in all validation cohorts. Thirdly, our study was based on data mining, and these findings should be further validated through biological experiments and clinical validations.
Despite these shortcomings, the RS5_G and RS_Lnc signatures constructed in this study performed well at predicting prognosis and grouping different prognostic groups.
In conclusion, an immune-related prognostic signature incorporating fewer genes was constructed in this study, which has comparable and even better at predicting and grouping abilities than those reported in other studies. Some immune-related lncRNAs are found to be correlated with the survival of NB patients, which are used to build the prognostic signature RS_Lnc, and it achieved good robustness. Collectively, this study not only provides prognostic signatures to help clinicians to choose the optimal therapeutic strategies but also sheds new light on research on NB treatment.

AUTHOR CONTRIBUTIONS
YL conceived of and directed the project. XZ designed the study, analyzed the data, and wrote the first draft of the manuscript. YT and HZ collected data and samples and prepared figures and tables. YZ and JC revised the manuscript critically for important intellectual content. LW reviewed the data and corrected the algorithm. All authors have read and approved the final draft for publication. All authors contributed to the article and approved the submitted version.

ACKNOWLEDGMENTS
We would like to acknowledge the Jilin Province Key Laboratory of Biometrics New Technology for its support of this project, and thanks to all researchers and patients who participated in the TARGET-NBL, R2, and GEO databases.