Machine learning methods revealed the roles of immune-metabolism related genes in immune infiltration,stemness,and prognosis of neuroblastoma

Abstract

BACKGROUND:

Immunometabolism plays an important role in neuroblastoma (NB). However, the mechanism of immune-metabolism related genes (IMRGs) in NB remains unclear. This study aimed to explore the effects of IMRGs on the prognosis, immune infiltration and stemness of patients with NB using machine learning methods.

METHODS:

R software (v4.2.1) was used to identify the differentially expressed IMRGs, and machine learning algorithm was used to screen the prognostic genes from IMRGs. Then we constructed a prognostic model and calculated the risk scores. The NB patients were grouped according to the prognosis scores. In addition, the genes most associated with the immune infiltration and stemness of NB were analyzed by weighted gene co-expression network analysis (WGCNA).

RESULTS:

There were 89 differentially expressed IMRGs between the MYCN amplification and the MYCN non-amplification group, among which CNR1, GNAI1, GLDC and ABCC4 were selected by machine learning algorithm to construct the prognosis model due to their better prediction effect. Both the K-M survival curve and the 5-year Receiver operating characteristic (ROC) curve indicated that the prognosis model could predict the prognosis of NB patients, and there was significant difference in immune infiltration between the two groups according to the median of risk score.

CONCLUSIONS:

We verified the effects of IMRGs on the prognosis, immune infiltration and stemness of NB. These findings could provide help for predicting prognosis and developing immunotherapy in NB.

Keywords

Neuroblastoma immune metabolism prognosis machine learning

1. Introduction

Neuroblastoma (NB) is a kind of neuroendocrine tumor derived from primitive cells of the sympathetic nervous system [1]. It is the most frequent tumor during infancy and accounts for about 15% of pediatric oncological deaths [2]. Despite spontaneous tumor regression observed in infants, its high mortality and malignancy remained a challenge for NB treatment. With the development of the novel discovery of the biological and genetic features of NB, various staging schemes for NB were used in the past. Age at diagnosis, stage, MYCN status, 1p and/or 11q status, INPC category, and individual histologic components were all considered to be an important risk classification factor [3]. With a prevalence rate of around 22% in NB patients [4], MYCN amplification is one of the alterations that correlates with poor prognosis and rapid tumor progression [5]. Although advanced progress has been made in NB treatment for low-risk and intermediate-risk patients, the five-year EFS and OS of high-risk cases remain low compared with other tumors, even after multimodal therapy, such as chemotherapy, surgery, radiotherapy, and Autologous Stem Cell Transplant (ASCT) [3, 6, 7]. Furthermore, late-onset adverse effects of current therapy have impacts on long-term survivors in all risk groups and are the barriers yet to be overcome [8, 9, 10]. Therefore, better treatments and further studies about tumorigenesis should be carried out to improve the prognosis of NB patients.

Immunotherapy is a revolutionary milestone in cancer treatment. With the development and thriving of modern immunology, researchers have established efficient and advanced therapies to bring benefits to cancer patients, including immune checkpoint inhibitors (ICIs) and adoptive cell transfer (ACT) [11]. However, different cancer may have different efficacies [12]. Pediatric tumors like NB seem to be insensitive to these implications. This result may be due to several facts, such as low mutation load and lack of MHC-1 expression, inactive tumor infiltration, and specific immune evasion strategies [13]. Although Jabbari et al. [14] demonstrated that the anti-GD2 antibody Dinutuximab is an effective NB treatment through molecular studies, the outcome is still frustrating. Nevertheless, the latest study once again shows us the hope of immunotherapy for NB. A phase I/II trial has revealed the promising activity of GD2-CART01. The GD2-CART01 cells expanded in 26 of 27 patients and persisted at detectable levels for a median of 3 months. Objective responses occurred in 17 patients (63%), including 9 (33%) with a complete response. Among patients treated at the recommended dose, 3-year event-free survival and overall survival were 36% and 60%, respectively [15]. As an emerging frontier, immunometabolism mainly focuses on the crosstalk between the immune cells and metabolism [16, 17]. Mathis et al. [18] suggested that immunometabolism play an important role in the field of obesity and type 2 diabetes initially. Recently, it has been discovered to be closely related to a variety of cancer [19, 20, 21]. It is generally acknowledged that the immune cells activate and infiltrate into tumors in response to their abnormal growth. To escape from immune surveillance and maintain sustainable growth, tumor cells compete with immune cells in nutrients, secrete metabolites or cytokines, and form tumor microenvironment (TME) after proliferation [22]. These distinct differences between cancer and normal cells require a series of reprogramming of cell metabolism. Warburg effect is one of the most classic, and it describes the elevated consumption of glucose in glycolysis rather than oxidative phosphorylation, even in the presence of physiological oxygen conditions in cancer cells [23]. Therefore, targeting the immune system to conquer the suppression caused by TME may be a vital opportunity for us to optimize the treatment of tumors. However, there are only studies about the impact of single immune or metabolic related genes on tumors at present [24, 25, 26]. Few noticed the importance of combined immune-metabolism related genes (IMRGs). Thus, a further understanding of immunometabolism may bring new insights and strategies for immunotherapy in NB.

Cancer stem cells (CSCs) are a group of tumor cells that keep generating new tissues consistently and lead to tumor growth [27]. With the revival of this concept, abundant evidence of the relationship between CSCs and tumor recurrence, metastasis or therapy resistance has been found [28, 29, 30, 31]. Moreover, it gives researchers innovative ideas for targets and therapies in the treatment of refractory tumors [32, 33]. A stemness index score (mRNAsi) is recently introduced to measure tumor stemness [34]. However, few studies concentrate on the correlation between CSCs and immunometabolism.

Machine Learning (ML), a major subfield of Artificial Intelligence (AI), enables computers to perform particular tasks after learning from a large amount of data automatically. Scientists have applied it to cancer research over the past decades. Models created by ML have contributed to the prediction of cancer types, prognosis, risk factors and treatment [35, 36]. In this study, we used datasets like GEO to establish a novel prognostic model based on IMRGs and analyze IMRG function in cancer stemness and immune infiltration. We will provide new insights into the role of IMRGs in NB, as well as the targets and clues for the diagnosis and treatment of NB.

2. Materials and methods

2.1 Data acquisition

“NB” and “Neuroblastoma” is used as a keyword to screen qualified gene data sets according to strict inclusion/exclusion criteria. The inclusion criteria are as follows: (1) Patients with sporadic NB; (2) clinical information of patients with NB, including MYCN status; (3) Test samples from the data set of human tissues. The exclusion criteria are: patients who have participated in clinical trials of drugs or other treatments. After filtering and comparison, according to strict inclusion/exclusion criteria, three gene datasets were selected from GEO (https://www.ncbi.nlm.nih.gov/GEO/): GSE120559 ( $n$ $=$ 208), GSE45547 ( $n=$ 649), GSE73517 ( $n=$ 105).MTAB-8248 ( $n=$ 223) and NB project dataset ( $n=$ 138) were selected from ArrayExpress and Target (https://ocg.cancer.gov/programs/target) databases. The data of Target and ArrayExpress (https://www.ebi.ac.uk/ biostudies/arrayexpress/) databases contain the overall survival rate and survival status of patients, which are used to establish and verify the prognosis model. We normalized the unstandardized data and corrected dataset GSE73517 that were not data-corrected using “limma” followed by log2 normalization. In addition, the R package “limma” was applied to combine the three datasets of GSE120559, GSE45547, and GSE73517 and remove the batch effect to obtain a new dataset ( $n=$ 932). 2483 immune-related genes were obtained from IMMPORT (https://www.immport.org/shared/home), and 325 metabolism-related genes were obtained from Gene Cards (https://www.genecards.org) website (Table S1). The genes with the correlation greater than 10 were selected using the keyword “metabolism”.

2.2 Differentially expressed analysis

The amplification of MYCN was one of the earliest biomarkers found in NB and is still considered to be one of the most reliable predictors of clinical manifestation [37]. As the amplification status of MYCN had a tremendous impact on the clinical phenotype of NB patients, the 932 NB samples after removal of batch effect were divided into two groups according to whether MYCN gene was amplified or not, including 173 MYCN-amplified patients and 759 MYCN non-amplified patients. The differentially expressed genes between the two groups were analyzed with the R package “limma” and the up-regulated and down-regulated genes were labeled. Performing intersection with immune-metabolism related genes obtain differentially expressed IMRGs.

2.3 Enrichment analysis

The R software “clusterProfiler” was used for Gene Ontology (GO) and Kyoto Encyclopedia of Genes and Genomics (KEGG) enrichment analysis of differentially expressed IMRGs. Additionally, we explored the quantitative relationship of genes among different enrichment pathways from website Metascape (https:// metascape.org/gp/#/main/step1).

2.4 Screening of prognostic genes

The “survival” package of R software was used for univariate Cox regression analysis of the differentially expressed IMRGs. STRING (version 11.0, https://string-db.org) database was utilized to figure out protein-protein interactions. A threshold of 0.15 was set, and the results were presented via Cytoscape software (v3.9.1). The “caret” package and the “randomForest” package were employed for Support Vector Machine Recursive Feature Elimination (SVM-RFE), and Random Forest (RF) algorithms to select the prognosis genes respectively. Support Vector Machine (SVM) is a kind of generalized linear classifier that binary classifies data according to supervised learning. The Recursive Feature Elimination (RFE) algorithm is used to screen the best gene to avoid overfitting. RF refers to a classifier that adopts multiple trees to train and predict samples. Ultimately, the screening results of prognostic genes were summarized by the “VennDiagram” package.

2.5 Identification of prognostic gene feature

Principal component analysis (PCA) was performed on the genotyping effect of prognostic genes in clinical phenotypes of patients with NB using R software “factoextra” and “FactoMineR” packages. We also analyzed the association of prognostic genes with immune checkpoints (ICs).

2.6 Construction of prognosis model

We chose the E-MTAB-8248 dataset as the train set ( $n=$ 223) and the dataset from the Target database as the test set ( $n=$ 137). The expression levels of the four selected prognostic genes were taken into account to predict the survival state and survival time of the patient. The model was constructed using the Least Absolute Shrinkage and Selection Operator (LASSO) algorithm in the package “glmnet” with 10,000 iterations. The risk score for each patient was calculated using the following formula:

$\textit{Risk Score}=\sum\textit{Coef}\left(i\right)\times\textit{Expr}\left(i\right)$

$\textit{Coef}\left(i\right)$ refers to coefficient representing each gene obtained by LASSO regression; $\textit{Expr}\left(i\right)$ refers to the expression level of the gene.

2.7 Validation of prognostic models

According to the median prognostic risk score acquired from prognostic model, the patients were divided into high-risk group and low-risk group. We used the “ROCR” and “timeROC” packages to calculate the 5 years’ Receiver operating characteristic (ROC) curve of the train set and the test set, and used the “survminer” package to calculate the K-M curve. Eventually, by means of the “rms” package, we built the nomogram and verified its accuracy.

2.8 Immune infiltration and stemness analysis

Stemness and immune infiltration are indicators of survival, classification and disease progression in cancer patients. Analysis of the two dimensions would help us to comprehend the mechanism of tumors’ action further. CIBERSORT is a method for characterizing cell composition of complex tissues from their gene expression profiles [38]. The patient’s immune infiltration was calculated by gene expression through the package “CIBERSORT”. The package “estimate” was applied to calculate the patient’s Stromal Score, Immune Score, ESTIMATE Score, Tumor Purity. The stemness index of a tumor can be predicted by gene expression. The one-class logistic regression machine learning algorithm (OCLR) extracted a set of transcriptomics and epigenetic features from untransformed pluripotent stem cells and their differentiated descendants to obtain a stemness index [34]. The model was built from the syn3156503 dataset in the “synapser” package, and the model was used to predict the stemness index of the train set and the test set, respectively.

2.9 Correlation analysis of IMRGs with immune score and stemness index

The operation of weighted gene co-expression network analysis (WGCNA) included clustering analysis of expression spectra and calculation of association between each cluster module and clinical phenotype [39]. We analyzed the correlation of the NB patient’s Stromal Score, Immune Score, ESTIMATE Score, Tumor Purity, and stemness index with the “WGCNA” package. The genes in the modules with most significant scores were analyzed for GO, KEGG, and Gene set enrichment analysis (GSEA) enrichment. By dint of Cytohubba plug-in in Cytoscape software [40], the genes with the most significant interactions were screened out.

2.10 Statistical analysis

Univariate Cox regression analysis, LASSO analysis, SVM-RFE and RF were selected to construct the prognosis model. The K-M survival curve was used to compare the prognosis differences between different risk groups. The ROC curve and Area Under Curve (AUC) values were used to assess the accuracy of the IMRGs-based prognostic model. All statistical analyses in this study were performed in R (v4.2.1). A $P$ -value of 0.05 was considered statistically significant.

3. Results

3.1 Study cohort

A total of 1323 patients were enrolled in this study from 5 data sets of 3 databases. After data cleaning, the clinical data of 1296 patients were left for statistical analysis. The Shapiro-Wilk normality test was employed to assess the normal distribution and levene test was used to assess the homogeneity of variance of clinical data [41]. Results showed that Age and INSS stage did not follow a normal distribution, and therefore, Wilcoxon and Kruskal-Wallis tests were used for their respective analyses. Conversely, the t-test and ANOVA were employed for other indices since they are normally distributed. The outcomes demonstrated that only the clinical data distribution of Risk group exhibited a significant difference (Table 1). The study was conducted in accordance with the Reporting Recommendations for Tumor Marker Prognostic Studies (REMARK) guidelines [42].

Table 1
Characteristics of study cohort

	MTAB-8248	GSE120559	GSE45547	GSE73517	NB project	Overall	$P$ value ${}^{\text{a}}$
	( $N=$ 223)	( $N=$ 208)	( $N=$ 649)	( $N=$ 105)	( $N=$ 111)	( $N=$ 1296)
Age
$<$ 1.5 years	104 (46.6%)	94 (45.2%)	414 (63.8%)	62 (59.0%)	26 (23.4%)	700 (54.0%)	0.818
$\geqslant$ 1.5 years	119 (53.4%)	114 (54.8%)	235 (36.2%)	43 (41.0%)	85 (76.6%)	596 (46.0%)
MYCN status
Amplified	46 (20.7%)	52 (25.1%)	93 (14.5%)	33 (31.4%)	$-$	224 (19.%)	0.153
Not amplified	176 (79.3%)	155 (74.9%)	550 (85.5%)	72 (68.6%)	$-$	953 (81.0%)
Unknown	1	1	6		111	113
INSS stage
1	29 (13.0%)	28 (13.5%)	153 (23.6%)	10 (9.5%)	(0%)	220 (17.0%)	0.296
2	39 (17.5%)	34 (16.3%)	113 (17.4%)	9 (8.6%)	1 (0.9%)	196 (15.1%)
3	36 (16.1%)	27 (12.9%)	91 (14.0%)	10 (9.5%)	5 (4.5%)	169 (13.0%)
4	89 (39.9%)	91 (43.8%)	214 (33.0%)	56 (53.3%)	85 (76.6%)	535 (41.3%)
4S	30 (13.5%)	28 (13.5%)	78 (12.0%)	20 (19.0%)	20 (18.0%)	176 (13.6%)
TERT rearrangement
TERT-RA	20 (9.0%)	21 (10.1%)	$-$	$-$	$-$	41 (9.5%)	0.069
Wild type	203 (91.0%)	187 (89.9%)	$-$	$-$	$-$	390 (90.5%)
Unkonwn			649	105	111	865
1p SCA
Normal	137 (67.2%)	$-$	$-$	57 (54.3%)	$-$	194 (62.8%)	0.302
SCA	67 (32.8%)	$-$	$-$	48 (45.7%)	$-$	115 (37.2%)
Unknown	19	208	649		111	987
Sex
Female	$-$	$-$	268 (43.8%)	$-$	45 (40.5%)	313 (43.3%)	0.655
Male	$-$	$-$	344 (56.2%)	$-$	66 (59.5%)	410 (56.7%)
Unknown	223	208	37	105		573
Risk group
High risk	$-$	$-$	$-$	56 (53.3%)	87 (78.4%)	143 (66.2%)	0.028
Intermediate risk	$-$	$-$	$-$	9 (8.6%)	10 (9.0%)	19 (8.8%)
Low risk	$-$	$-$	$-$	40 (38.1%)	14 (12.6%)	54 (25.0%)
Unkown	223	208	649			1080

Abbreviations: TERT-RA, TERT rearrangement; SCA, segmental chromosome aberrations.

Figure 1.

Identification of differentially expressed genes. (A) Data-uncorrected expression data. (B) Data-corrected expression data. C. Volcanic map of differentially expressed genes. D. Heatmap of gene expression for 932 NB patients.

3.2 Identification of IMRGs

To identify differentially expressed IMRGs, we analyzed three datasets, GSE120559 ( $n=$ 208), GSE45547 ( $n=$ 649) and GSE73517 ( $n=$ 105). Since the datasets were presented with uncalibrated data (Fig. 1A), we calibrated the data using R software (Fig. 1B). The calibrated data were combined, and the batch effect was removed by virtue of the function “RemoveBatchEffect” of the R package “limma”. Next, we obtained a new expression data sample with a size of 932 after adjusting the abnormal expression data. This new sample was divided into two groups according to whether MYCN was amplified, the MYCN amplification group ( $n=$ 173) and the MYCN non-amplification group ( $n=$ 759). As for differentially expressed genes, a total of 811 differentially expressed genes were identified using the “limma” package, of which 464 differentially up-regulated genes and 346 differentially down-regulated genes (Fig. 1C–D). There were 89 differentially expressed IMRGs, of which 67 were up-regulated and 22 were down-regulated.

Figure 2.

Enrichment analysis of differentially expressed IMRGs. (A) GO enrichment analysis. The outermost ring is the enrichment pathway; The second circle is the number of genes enriched into the pathway; The third circle is the gene expression; The innermost circle is the enrichment significant degree. (B) KEGG enrichment analysis. On the right are the enrichment pathways, and on the left are the top 15 most significant genes involved in enriching each pathway. (C) GO enrichment BP pathway results. GO Cluster generates a circular dendrogram of the clustering using by default euclidean distance and average linkage. The inner ring displays the color coded logFC while the outside one encodes the assigned terms to each gene. (D) GO enrichment results of Metascape database. The points represent the enrichment pathway. The larger the point, the more genes enriched into the pathway. The larger the line between the points, the more common genes between the points.

3.3 Enrichment analysis

To further explore the biological effects of the differentially expressed IMRGs, enrichment analysis was performed on the 89 genes (Fig. 2A). GO enrichment analysis showed that in biological processes (BP), the differential genes were enriched in 1116 pathways, including lymphocyte differentiation, T cell differentiation, and monocyte differentiation (Fig. 2B); In cellular components (CC), differential genes were enriched in 10 pathways, including the alpha-beta cell receptor complex; In molecular functions (MF), the differential genes were enriched in 129 pathways, including the cytokine activity, Signaling Receptor Activity. The KEGG enrichment analysis manifested differentially expressed IMRGs enrichment in hematopoietic cell lines, chemokine signaling pathways, cytokine interactions, and T cell receptors (Fig. 2C). The analysis from website Metascape showed that leukocyte activity had the largest number of differentially expressed IMRGs and relevant pathways (Fig. 2D).

3.4 Screening of prognostic genes

Figure 3.

Screening of prognostic genes. (A) The optimal number of parameters using SVM-RFE algorithm. (B) The optimal value using random forest algorithm. (C) The prognostic genes using random forest algorithm. (D) Top 20 genes most related to prognosis in univariate Cox regression analysis results. (E) Protein-protein interaction network diagram of top 15 and other 74 genes in univariate Cox regression analysis results. (F) Prognosis genes obtained by combining the SVM-RFE algorithm, the RF algorithm and the univariate Cox regression.

In order to screen out the genes most related to the prognosis of patients with NB to construct a prognosis model, we adopted machine learning method to select genes. We use the function “rfe” in the “caret” package of R software, and set the train method as “svmRadial”. The results showed that the expression data of seven genes selected best predicted the prognosis of the patients (Fig. 3A). Using the “randomForest” package of R software, the optimal number of decision trees and parameters were calculated (Fig. 3B), and these optimal parameters were then used to construct a model that obtained 37 genes (Fig. 3C). Furthermore, 20 genes most associated with the prognosis of patients with NB were selected using univariate Cox regression (Fig. 3D), and the protein-protein interaction of these genes with the differentially expressed IMRGs was shown in a network diagram (Fig. 3E). Finally, we screened out four prognostic genes by intersecting the genes selected by SVM-RFE, RandomFrorest and univariate Cox regression, of which CNR1 and GLDC were metabolism-related genes, while GNAI1 and ABCC4 were immune-related genes (Fig. 3F).

Figure 4.

Correlation analysis between prognostic genes and clinical features and immune checkpoints of NB patients. (A–C) Clustering effect of four prognostic genes on TERT gene mutation, MYCN gene amplification and tumor staging using PCA. (D) Correlation of four prognostic genes with five immune checkpoints in the train set. The color of the line represents the correlation between the two, ranging from $-$ 1 to 1. (E) Correlation of four prognostic genes with five immune checkpoints in the test set. (F–I) K-M curves of four prognostic genes and prognosis of NB patients. The top graph is K-M curve, the number in the middle table is the number of surviving patients, and the line in the bottom table is that patients have died or withdrawn from the study.

To further explore the role of the four prognostic genes in patients with NB, PCA was used to identify the effect of classification of 4 prognostic genes on the clinical features of NB patients (Fig. 4A–C). The correlation between immune checkpoints (ICs) CD27, CD40, CD48, CTLA4, and LAG3 and the expression of prognostic genes showed that CNR1, GNAI1 tended to be negatively correlated with ICs’ expression, while GLDC and ABCC4 tended to be positively correlated with ICs’ expression in train set (Fig. 4D–E). The Kaplan-Meier (K-M) survival curves of 4 prognostic genes and NB patients were plotted using R software (Fig. 4F–I). The results showed that high expression of GLDC and ABCC4 might lead to poor prognosis in NB patients, while low expression of CNR1 and GNAI 1 might lead to poor prognosis in NB patients.

3.5 Construction and verification of prognostic model

Figure 5.

Establishment of prognostic model. (A) Process of establishing prognostic model using LASSO regression. Each curve in the graph represents the change trajectory of the coefficient of each independent variable, with the ordinate being the value of the coefficient, the lower abscissa being log( $\lambda)$ , and the upper abscissa being the number of non-zero coefficients in the model at this time. With the change of $\lambda$ value, the later the coefficient is compressed to 0, the more important it is. (B) There are two numerical dashed lines in the figure, the one with the lowest error on the left and the one with fewer features on the right. (C, D) Distribution of risk scores of NB patients with different survival states in the train set and the test set. (E, F) K-M curves between train set and NB patients with higher risk score and lower risk score in the test set. (G, H) 5-year ROC curves of prognostic models in the train set and the test set.

Figure 6.

Analysis of prognostic risk scores and clinical features. (A) Correlation between clinical information in train set and test set of univariate Cox regression analysis, risk score and prognosis of NB patients. (B) Nomograms constructed using clinical information and risk scores. (C) Accuracy of nomograms in predicting the prognosis of patients with NB. The abscissa is the Predicted risk, and the ordinate is the Observed risk, ranging from 0 to 1, which can be understood as the incident rate (percentage). The dashed diagonal line is the reference line, that is, the predicted value = the actual value. The red line is a curve fitting line, and the colored parts on both sides are 95%CI. (D, E) Mulberry plots show the distribution of clinical information and risk scores in the train and test sets.

Aiming to further analyze the prognostic role of prognostic genes in NB patients, a prognostic prediction model was constructed using 4 prognostic genes and LASSO regression. The E-MTAB-8248 data set was selected as the train set ( $n=$ 223) and the data set from the Target database as the test set ( $n=$ 137) and the optimal $\lambda$ value was selected for prognostic model (Fig. 5A–B). Prognostic model parameter information was obtained after 10,000 iterations (Table S2). The model can be used to calculate the prognostic risk scores of each NB patient. It can be found that the prognostic scores of NB patients with different survival outcomes showed significant differences in the train set and the test set (Fig. 5C–D). The median prognostic risk scores obtained using the prognostic model divided the patients into high-risk and low-risk groups. The K-M survival curves of different groups also showed significant differences in the train set and the test set (Fig. 5E–F). The 5-year ROC curve could estimate the accuracy of the prognosis model in predicting the prognosis of NB patients within 5 years. The AUC of the 5-year ROC curves of the train set were 0.706, 0.771, 0.776, 0.810, and 0.825 (Fig. 5G), respectively. The AUC of the five-year ROC curves in the test set were 0.565, 0.600, 0.662, 0.646, and 0.650, respectively (Fig. 6H). Univariate Cox regression was used to analyze the prognostic factors of patients. The risk score was a significant prognostic factor in the train set but not in the test set (Fig. 6A). To facilitate the prognosis evaluation of patients, nomograms were constructed using clinical information and prognostic risk scores (Fig. 6B) with good predictive accuracy (Fig. 6C). Finally, we elucidated the correlation between survival information for NB patients, age, MYCN amplification status, tumor stage, TERT mutation status, sex, ethnicity, and risk score (Fig. 6D–E).

Figure 7.

Analysis of immune infiltration. (A, B) Percentage of immune cells in train set and test set. (C, D) Difference of immune cell infiltration between NB patients with higher risk score and lower risk score in train set and test set. (E, F) Stromal Score, Immune Score, ESTIMATE Score, Tumor Purity in the train set and test set.

Figure 8.

Correlation analysis between prognosis gene and prognosis of NB patients and immune infiltration. (A) Mantel’test analysis of prognosis of NB patients and immune infiltration. (B, C) Risk scores in the train and test sets were higher than the stemness index for NB patients with low risk scores. (D, E) Radar diagram of correlation between four prognostic genes and immune score, matrix score, estimated score and stemness index in the train set and test set. (F) Correlation network diagram of four prognostic genes and 22 immune cells.

3.6 Analysis of immune infiltration and stemness of IMRGs

The analysis of immune infiltration and stemness in patients with NB can provide new ideas for the pathogenesis and treatment of NB. The immune infiltration of NB patients was calculated using the R package “CIBERSORT” (Fig. 7A–B). By grouping the prognosis models, we compared the differences in immune infiltration between the two groups, which shows that in the train set, the infiltration of Plasma cells, T.cells.CD4.memory.resting, T.cells.CD4.memory.activated, T.cells.follicular.helper and Neutrophils was significantly different between the two groups, while B.cells.naive, B.cells.memory and T.cells.regulatory and Tregs differ significantly in the test set (Fig. 7C–D). The “estimate” package was used to calculate the Stromal Score, Immune Score,ESTIMATE Score and Tumor Purity for patients with NB. The relationship of these scores to the clinical features of NB patients was also shown (Fig. 7E–F). The correlation between immune infiltration cells and the patient’s prognosis was analyzed via mantel’test, which showed that T.cells.CD4.naïve was most relevant to the prognosis of patients with NB (Fig. 8A). Besides, the stemness index of NB patients was calculated by the stemness index calculation model (Fig. 8B–C), which revealed a trend towards higher stemness index in NB patient with a high risk score. At last, we analyzed the association of the 4 prognostic genes with the Stromal Score, Immune Score, ESTIMATE Score, Tumor Purity and stemness index (Fig. 8D–F).

3.7 Screening of key genes for immune infiltration and stemness

Figure 9.

Genes related to Stromal Score, Immune Score, ESTIMATE Score, Tumor Purity and stemness index screened using WGCNA. (A) Determination of soft threshold and power values. To determine the optimal soft threshold, we mainly refer to the left figure, that is, the scale-free fit index (Y axis) under different soft thresholds (X axis). The red line indicates the value of the subjectively selected scale-free fitting index, as shown in figure 0.9. When the scale-free fitting index is 0.9, the minimum soft threshold for constructing scale-free network is 4, so 4 can be selected as the optimal soft threshold for subsequent analysis. The figure on the right shows the network connectivity under different soft thresholds. (B) Cluster tree diagram of the patient’s gene. Each branch in the diagram represents a gene, and each color below represents a co-expression module. (C) clustering tree diagram of NB patients and correlation heat diagram of clinical characteristics. (D) Heatmap of modules correlation with Stromal Score, Immune Score, ESTIMATE Score, Tumor Purity and stemness index. (E) GO enrichment analysis of Brown and Turquoise modules.

Figure 10.

Enrichment analysis of the genes selected by WGCNA. (A) KEGG Analysis Results. (B) GSEA Results. (C) Protein-protein interaction network diagram for the 10 genes with the highest MCC scores calculated by the “Cytohubba” plug-in.

The identification of key genes may provide potential targets for the treatment of patients with NB. The IMRGs were clustered into different modules, each containing multiple genes, using the R software “WGCNA” package. To begin with, the best soft threshold was found and set to 5 (Fig. 9A) and the adjacency matrix and TOM are established accordingly (Fig. 9B). Subsequently, the clustering of NB patients and the correlation analysis between immune infiltration and stemness were performed (Fig. 9C). Then, a total of 9 modules were identified. Correlation Analysis between each module and the immune infiltration and stemness showed that the module brown and turquoise had a significant positive correlation with the Stromal Score, Immune Score, ESTIMATE Score, and a significant negative correlation with the Tumor Purity and stemness index (Fig. 9D). We selected genes from the brown module and turquoise module, and conducted gene enrichment analysis using GO enrichment analysis, KEGG enrichment analysis, and GSEA. GO enrichment results showed that the two modular genes were mainly enriched in the cytokine-mediated signaling pathway and receptor ligand activity (Fig. 9E). The KEGG enrichment results showed that the two modular genes were mainly enriched in the cytosine receptor interaction, PI3K-Akt signaling pathway, Chemokine signaling pathway (Fig. 10A). GSEA results showed genes enriched in lymphocyte-activated, lymphocyte-mediated immunity (Fig. 10B). The Cytoscape software Cytohubba plug-in was used to select the key genes. According to results analyzed by the ‘Cytohubba’, the top 10 genes with the highest MCC scores were selected as the hub genes (Fig. 10C).

4. Discussion

NB is a cancer of the sympathetic system and represents the most common extracranial solid tumor in childhood [43, 44]. Although spontaneous regression is observed in certain cases, others with MYCN amplification metastasize rapidly and result in a large number of deaths [45]. These high-risk metastatic NBs have characteristics of cold tumors since the native immune response is considered inactive in the anti-tumor process. With a better understanding of NB, novel immunotherapy using anti-GD2 antibody has been invented, opening a new era for NB treatment [46]. The survival rate of high-risk NBs has increased to 50% while most patients suffer from adverse effects after multimodal treatment like chemotherapy, surgical resection or radiotherapy [47]. Moreover, relapse and substantial neuropathic pain caused by anti-GD2 antibody are still challenges [48]. Cancer metabolism was a concept popularized by Otto Warburg nearly 100 years ago. Nowadays, researchers realize that it is not only comprised of intrinsic metabolic alterations of cancer cells but involved with metabolic demands of immune and other normal cells in the tumor microenvironment (TME) [49]. The interaction of metabolic reprogramming between cancer and immune cells is regarded as a pivotal determinant in anti-tumor immune response [50]. Therefore, a meticulous study on IMRGs in NB may give us an opportunity to have insights into tumorigenesis and help us develop innovative immunotherapy and improve the prognosis of NB patients.

In this study, we first identified differentially expressed IMRGs and performed enrichment analysis on them. Then, machine learning algorithm was used to screen out the possible prognostic genes and LASSO regression was used to establish a prognostic model, and the prognostic risk score was calculated using the prognostic model. In addition, as much progress has been made in tumor studies with tumor stemness and immune infiltration, the genes most significantly associated with NB stemness index and immune score were identified with WGCNA.

Enrichment analysis results showed that the differentially expressed IMRGs were enriched in the lymphocyte differentiation, T cell differentiation, and monocyte differentiation pathways. Significant inhibition of CD4+T lymphocyte differentiation has been shown in 9,10-secosteroids [51]; dendritic cell 2 analog iloprost treated with prostaglandin I promotes antigen-specific regulatory T-cell differentiation in mice [52]; the QKI subtype controls cell differentiation by regulating alternative splicing, mRNA stability, and translation [53]. What’s more, analysis of immune infiltration also revealed significant differences in CD4+T cell infiltration between the two groups with different risk scores. The above results indicated that there were specific drugs to regulate immune cell differentiation, which might also be a potential target for drug treatment of NB.

We constructed a prognostic model based on four IMRGs (CNR1, GNAI1, GLDC, and ABCC4). Results illustrated through the KM survival curve show that low CNR1 and GNAI1 expression or high expression of GLDC and ABCC4 is related to poor prognosis in NB patients. CB1, expressed by CNR1, is a cannabinoid receptor in the endocannabinoid system (ECS) and a member of the G-protein coupled receptor (GPCR) family [54]. Researchers found that CB1 receptors are expressed by various immune cells and may regulate liver metabolism [55]. Its ligand endocannabinoid (EC) impacts tumorigenesis and immune response [56]. Moreover, CB1-specific agonists suppress proliferation and vascularization while enhancing apoptosis in non-small cell lung cancer (NSCLC) by inhibiting phosphorylation of AKT [57], indicating CB1 may be a novel target in NSCLC treatment. A similar correlation between CB1 expression and the AKT signaling pathway is also observed in prostate cancer [58]. However, no study has elucidated the mechanism of CNR1 in NB development. Further study should be investigated on the relationship between CNR1 and NB. The guanine nucleotide-binding protein G(i), $\alpha$ -1 subunit (GNAI1), is one of the signaling components of G proteins [59]. Interestingly, it plays an important role in VEGF-induced angiogenesis and upregulates in patients with PDR [60]. Along with GNAI2 and GNAI3, GNAI1 is highly expressed in immune cells and restrains colon tumor development in mice by inhibiting the IL-6-mediated pathway in CD4+/CD11b+ DCs [61]. Furthermore, the downregulation of GNAI1 allows migration and invasion of HCC cells [59]. Therefore, GNAI1 may be a protective factor in tumorigenesis. GLDC is a metabolic oncogene confirmed recently. As the product of GLDC expression, glycine decarboxylase (GLDC) is part of a multienzyme complex (the glycine cleavage system, GCS) and catalyzes the process of glycine decomposition [62]. Glycine metabolism is an innovative target in hematologic malignancies [63, 64]. It has been discovered that a lack of GLDC causes neural tube defects (NTDs) by decreasing the one-carbon unit supply in mitochondrial folate metabolism [65]. Studies have found that GLDC promotes tumor growth in several cancers, such as HCC, NSCLC, and lung cancer [66, 67, 68]. Furthermore, the expression of GLDC increases significantly and is identified as a transcriptional target crucial for cancer cell proliferation in MYCN-amplified NB [69]. The ATP-binding cassette (ABC) superfamily consists of a series of transmembrane proteins and is acknowledged to be responsible for multi-drug resistance through drug efflux in many tumors [70, 71]. ABCC4, also known as MRP4, is the shortest one of the ABCC family and is ubiquitous in tissue cells [72]. Knockdown of ABCC4 has been found to result in suppression of proliferation and migration in different tumor cells [73, 74, 75, 76]. Moreover, a study has shown that high ABCC4 expression is independently correlated with poor prognosis in NB, other than MYCN amplification, age at diagnosis, or disease stage, elucidating that it could be a prognostic biomarker [77]. Although deletion of ABCC4 improves the sensitivity to irinotecan in NB patients and slows tumor growth, it does not decrease tumor incidence fundamentally [78]. Whereas, related research has found that ABCC4 and ABCC1 alteration influences tumor progression even when no chemotherapy treatment is given. This discovery indicates that the impacts of ABCC4 on NB are mediated through other pathways except for cytotoxic drug efflux [79]. In conclusion, these results may contribute to NB treatment and serve as potential targets in the future.

Cancer stem cells (CSCs) can maintain tumors by manipulating genetic and non-genetic factors to metastasize, resist treatment, and maintain the tumor microenvironment. Understanding the key features and mechanisms of CSCs is conducive to improve patients’ prognosis by improving prognostic models and treatments [80]. In recent years, there have been many reports on the interaction between immunity and stemness in tumors. For example, cancer cells evade the immune system by epigenetic reprogramming to increase stemness [81]. In our study, WGCNA was used to identify the gene modules with the most significant correlation with NB stemness index and immune score, and the selected gene modules were found to be significantly positively correlated with Stromal Score, Immune Score, and significantly negatively correlated with stemness index and Tumor Purity, indicating that the genes contained in these modules might lead to a good prognosis for patients.

5. Conclusion

In conclusion, we developed an IMRGs-based model for the prognosis of NB patients and has a brilliant effect of prediction. Further, the pivotal genes which also play a crucial role in the immunity and stemness of NB were screened out, which may be used as new sites for the research on the mechanism of NB in the future. The above results illustrate the primary role of IMRGs in NB. However, these conclusions are based on the analysis of data from public databases, and further large-scale experiments are needed to verify them in further studies.

Abbreviations

NB:	Neuroblastoma
IMRG:	Immune-metabolism related gene
GEO:	Gene Expression Omnibus
DEG:	Differentially expressed gene
KEGG:	Kyoto Encyclopedia of Genes and Genomes
CIBERSORT:	Cell type identified by estimating relative subsets of RNA transcripts
GSEA:	Gene Set Enrichment Analysis
LASSO:	Least absolute shrinkage and selection operator
KM:	Kaplan–Meier
ROC:	Receiver operating characteristic
AUC:	Area under the curve
OS:	Overall survival
PCA:	Principal component analysis
DNA:	Deoxyribonucleic acid

Data and code availability

This study analyzed publicly available data sets. 3 datasets were selected from GEO (https://www.ncbi.nlm.nih.gov/GEO/): GSE120559 ( $n=$ 208), GSE45547 ( $n=$ 649), GSE73517 ( $n=$ 105). MTAB-8248 ( $n=$ 223) and NB project dataset ( $n=$ 138) was selected from ArrayExpress and Target (https://ocg.cancer.gov/programs/target) databases. 2483 immune-related genes were obtained from IMMPORT (https://www.immport.org/shared/home), and 325 metabolism-related genes were obtained from Gene Cards (https://www.genecards.org) website. The codes used in our study are available on reasonable request from the author of the communication. The code to reproduce the analyses in this manuscript using our R package are also available via GitHub (https://github.com/Jianhua-Mu/The_main_codes).

Author contributions

J.M. contributed to research design and data analysis. J.G. assisted editing and sorted data. P.L. and M.Z. assisted in the revision of manuscripts and documentation of references. K.W directed the project, revised the manuscript, and provided writing guidance. All authors contributed to the revision of the manuscript and read and approved the submitted version.

Funding

This study was supported by the Guangdong Basic and Applied Basic Research Foundation of China (Grant No. 2019 A1515011086) and College Students’ Entrepreneurship and Innovation Training Program (Grant No. 2022 12121229).

Ethics approval and consent to participate

Not applicable.

Consent for publication

Not applicable.

Competing interests

The authors declare that they have no competing interests.

Supplementary data

The supplementary files are available to download from https://dx-doi-org.web.bisu.edu.cn/10.3233/CBM-230119.

sj-xlsx-1-cbm-10.3233_CBM-230119.xlsx - Supplemental material

Supplemental material, sj-xlsx-1-cbm-10.3233_CBM-230119.xlsx

sj-xlsx-2-cbm-10.3233_CBM-230119.xlsx - Supplemental material

Supplemental material, sj-xlsx-2-cbm-10.3233_CBM-230119.xlsx

Footnotes

Acknowledgments

Not applicable.

References

Brodeur

G.M.

, Neuroblastoma: biological insights into a clinical enigma, Nat Rev Cancer 3(3) (2003), 203–216.

Maris

J.M.

Hogarty

M.D.

Bagatell

and Cohn

S.L.

, Neuroblastoma, Lancet 369(9579) (2007), 2106–2120.

Irwin

M.S.

et al., Revised Neuroblastoma Risk Classification System: A Report From the Children’s Oncology Group, J Clin Oncol 39(29) (2021), 3229–3241.

Brodeur

G.M.

Maris

J.M.

Yamashiro

D.J.

Hogarty

M.D.

and White

P.S.

, Biology and genetics of human neuroblastomas, J Pediatr Hematol Oncol 19(2) (1997), 93–101.

Seege

R.C.

et al., Association of multiple copies of the N-myc oncogene with rapid progression of neuroblastomas, N Engl J Med 313(18) (1985), 1111–1116.

Ladenstein

et al., 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(4) (2017), 500–514.

Park

J.R.

et al., Effect of Tandem Autologous Stem Cell Transplant vs Single Transplant on Event-Free Survival in Patients With High-Risk Neuroblastoma: A Randomized Clinical Trial, Jama 322(8) (2019), 746–755.

Friedman

D.N.

and Henderson

T.O.

, Late Effects and Survivorship Issues in Patients with Neuroblastoma, Children (Basel) 5(8) (2018), 102.

Cohen

L.E.

et al., Late effects in children treated with intensive multimodal therapy for high-risk neuroblastoma: high incidence of endocrine and growth problems, Bone Marrow Transplant 49(4) (2014), 502–508.

10.

Gurney

J.G.

Tersak

J.M.

Ness

K.K.

Landier

Matthay

K.K.

and Schmidt

M.L.

, Hearing loss, quality of life, and academic problems in long-term neuroblastoma survivors: a report from the Children’s Oncology Group, Pediatrics 120(5) (2007), e1229–1236.

11.

Pan

et al., Next-generation immuno-oncology agents: current momentum shifts in cancer immunotherapy, J Hematol Oncol 13(1) (2020), 29.

12.

Zhang

and Zhang

, The history and advances in cancer immunotherapy: understanding the characteristics of tumor-infiltrating immune cells and their therapeutic implications, Cell Mol Immunol 17(8) (2020), 807–821.

13.

Wienke

Dierselhuis

M.P.

Tytgat

G.A.M.

Künkele

Nierkens

and Molenaar

J.J.

, The immune landscape of neuroblastoma: Challenges and opportunities for novel therapeutic strategies in pediatric oncology, Eur J Cancer 144 (2021), 123–150.

14.

Jabbari

Hanaei

and Rezaei

, State of the art in immunotherapy of neuroblastoma, Immunotherapy 11(9) (2019), 831–850.

15.

Killock

, CAR T cells induce durable remission of neuroblastoma, Nat Rev Clin Oncol 20(6) (2023), 354.

16.

O’Neill

L.A.

Kishton

R.J.

and Rathmell

, A guide to immunometabolism for immunologists, Nat Rev Immunol 16(9) (2016), 553–565.

17.

Chen

Wang

Tang

and Kang

, ACOD1 in immunometabolism and disease, Cell Mol Immunol 17(8) (2020), 822–833.

18.

Mathis

and Shoelson

S.E.

, Immunometabolism: an emerging frontier, Nat Rev Immunol 11(2) (2011), 81.

19.

Zhang

Lou

Bai

X.L.

and Liang

T.B.

, Immuno- metabolism: A novel perspective of liver cancer microenvironment and its influence on tumor progression, World J Gastroenterol 24(31) (2018), 3500–3512.

20.

Yang

Chen

Gao

and Wang

, Tumor Immuno- metabolism Characterization in Ovarian Cancer With Prognostic and Therapeutic Implications, Front Oncol 11 (2021), 622752.

21.

Zhu

et al., Clinical Significance and Immunometabolism Landscapes of a Novel Recurrence-Associated Lipid Meta- bolism Signature In Early-Stage Lung Adenocarcinoma: A Comprehensive Analysis, Front Immunol 13 (2022), 783495.

22.

Traba

Sack

M.N.

Waldmann

T.A.

and Anton

O.M.

, Immunometabolism at the Nexus of Cancer Therapeutic Efficacy and Resistance, Front Immunol 12 (2021), 657293.

23.

Koppenol

W.H.

Bounds

P.L.

and Dang

C.V.

, Otto Warburg’s contributions to current concepts of cancer metabolism, Nat Rev Cancer 11(5) (2011), 325–337.

24.

et al., Identification of Iron Metabolism-Related Gene Signatures for Predicting the Prognosis of Patients With Sarcomas, Front Oncol 10 (2020), 599816.

25.

Shang

Huang

and Guo

, Identification of an Metabolic Related Risk Signature Predicts Prognosis in Cervical Cancer and Correlates With Immune Infiltration, Front Cell Dev Biol 9 (2021), 677831.

26.

Chen

Z.Y.

Zhou

G.Q.

and Sun

, An Immune-Related Gene Prognostic Index for Head and Neck Squamous Cell Carcinoma, Clin Cancer Res 27(1) (2021), 330–341.

27.

Batlle

and Clevers

, Cancer stem cells revisited, Nat Med 23(10) (2017), 1124–1134.

28.

Huang

et al., Stem cell programs in cancer initiation, progression, and therapy resistance, Theranostics 10(19) (2020), 8721–8743.

29.

Wang

et al., Single-Cell Analyses Reveal Mechanisms of Cancer Stem Cell Maintenance and Epithelial-Mesenchymal Transition in Recurrent Bladder Cancer, Clin Cancer Res 27(22) (2021), 6265–6278.

30.

Bai

Beretov

Graham

and Li

, Cancer stem cell in breast cancer therapeutic resistance, Cancer Treat Rev 69 (2018), 152–163.

31.

Xie

X.P.

et al., Quiescent human glioblastoma cancer stem cells drive tumor initiation, expansion, and recurrence following chemotherapy, Dev Cell 57(1) (2022), 32–46e38..

32.

Quaglino

Conti

and Cavallo

, Breast cancer stem cell antigens as targets for immunotherapy, Semin Immunol 47 (2020), 101386.

33.

Han

et al., Cancer stem cell-targeted bio-imaging and chemotherapeutic perspective, Chem Soc Rev 49(22) (2020), 7856–7878.

34.

Malta

T.M.

et al., Machine Learning Identifies Stemness Features Associated with Oncogenic Dedifferentiation, Cell 173(2) (2018), 338–354e315.

35.

Le Glaz

et al., Machine Learning and Natural Language Processing in Mental Health: Systematic Review, J Med Internet Res 23(5) (2021), e15708.

36.

Kourou

Exarchos

T.P.

Exarchos

K.P.

Karamouzis

M.V.

and Fotiadis

D.I.

, Machine learning applications in cancer prognosis and prediction, Comput Struct Biotechnol J 13 (2015), 8–17.

37.

Park

J.R.

et al., Children’s Oncology Group’s 2013 blueprint for research: neuroblastoma, Pediatr Blood Cancer 60(6) (2013), 985–993.

38.

Newman

A.M.

et al., Robust enumeration of cell subsets from tissue expression profiles, Nat Methods 12(5) (2015), 453–457.

39.

Langfelder

and Horvath

, WGCNA: an R package for weighted correlation network analysis, BMC Bioinformatics 9 (2008), 559.

40.

Shannon

et al., Cytoscape: a software environment for integrated models of biomolecular interaction networks, Genome Res 13(11) (2003), 2498–2504.

41.

Shapiro

S.S.

and Wilk

M.B.

, An Analysis of Variance Test for Normality (Complete Samples), Biometrika 52(3/4) (1965), 591–611.

42.

McShane

L.M.

Altman

D.G.

Sauerbrei

Taube

S.E.

Gion

and Clark

G.M.

, REporting recommendations for tumour MARKer prognostic studies (REMARK), Br J Cancer 93(4) (2005), 387–391.

43.

Seeger

R.C.

, Immunology and immunotherapy of neuroblastoma, Semin Cancer Biol 21(4) (2011), 229–237.

44.

A.L.

et al., Anti-GD2 antibody with GM-CSF, interleukin-2, and isotretinoin for neuroblastoma, N Engl J Med 363(14) (2010), 1324–1334.

45.

Ambros

P.F.

et al., International consensus for neuroblastoma molecular diagnostics: report from the International Neuroblastoma Risk Group (INRG) Biology Committee, Br J Cancer 100(9) (2009), 1471–1482.

46.

Cheung

N.K.

and Dyer

M.A.

, Neuroblastoma: developmental biology, cancer genomics and immunotherapy, Nat Rev Cancer 13(6) (2013), 397–411.

47.

Park

J.A.

and Cheung

N.V.

, Targets and Antibody Formats for Immunotherapy of Neuroblastoma, J Clin Oncol 38(16) (2020), 1836–1848.

48.

Anderson

Majzner

R.G.

and Sondel

P.M.

, Immunotherapy of Neuroblastoma: Facts and Hopes, Clin Cancer Res 28(15) (2022), 3196–3206.

49.

Dang

C.V.

, Cancer Metabolism: the Known, Unknowns, Biochim Biophys Acta Rev Cancer 1870(1) (2018), 1.

50.

Xia

et al., The cancer metabolic reprogramming and immune response, Mol Cancer 20(1) (2021), 28.

51.

Sun

Y.L.

Tang

Zheng

G.L.

and Zhang

, Immunosuppressive 9,10-Secosteroids from the Gorgonian Verrucella umbraculum Collected in the South China Sea, J Nat Prod 84(5) (2021), 1671–1675.

52.

Wong

T.H.

et al., Dendritic cells treated with a prostaglandin I(2) analog, iloprost, promote antigen-specific regulatory T cell differentiation in mice, Int Immunopharmacol 79 (2020), 106106.

53.

Neumann

D.P.

Goodall

G.J.

and Gregory

P.A.

, The Quaking RNA-binding proteins as regulators of cell differentiation, Wiley Interdiscip Rev RNA 13(6) (2022), e1724.

54.

Braile

Marcella

Marone

Galdiero

M.R.

Varricchi

and Loffredo

, The Interplay between the Immune and the Endocannabinoid Systems in Cancer, Cells 10(6) (2021), 1282.

55.

Liu

L.Y.

et al., Cannabinoid receptor signaling regulates liver development and metabolism, Development 143(4) (2016), 609–622.

56.

Śledziński

Zeyland

Słomski

and Nowak

, The current state and future perspectives of cannabinoids in cancer biology, Cancer Med 7(3) (2018), 765–775.

57.

Preet

et al., Cannabinoid receptors, CB1 and CB2, as novel targets for inhibition of non-small cell lung cancer growth and metastasis, Cancer Prev Res (Phila) 4(1) (2011), 65–75.

58.

Cipriano

Häggström

Hammarsten

and Fowler

C.J.

, Association between cannabinoid CB1 receptor expression and Akt signalling in prostate cancer, PLoS One 8(6) (2013), e65798.

59.

Yao

et al., GNAI1 Suppresses Tumor Cell Migration and Invasion and is Post-Transcriptionally Regulated by Mir-320a/c/d in Hepatocellular Carcinoma, Cancer Biol Med 9(4) (2012), 234–241.

60.

Sun

et al., Gαi1 and Gαi3mediate VEGF-induced VEGFR2 endocytosis, signaling and angiogenesis, Theranostics 8(17) (2018), 4695–4709.

61.

Z.W.

et al., GNAI1 and GNAI3 Reduce Colitis-Associated Tumorigenesis in Mice by Blocking IL6 Signaling and Down-regulating Expression of GNAI2, Gastroenterology 156(8) (2019), 2297–2312.

62.

M.K.

Zhang

W.C.

Lim

and Yew

W.S.

, Glycine decarboxylase is an unusual amino acid decarboxylase involved in tumorigenesis, Biochemistry 53(5) (2014), 947–956.

63.

Liu

et al., Loss of PRMT7 reprograms glycine metabolism to selectively eradicate leukemia stem cells in CML, Cell Metab 34(6) (2022), 818–835e817..

64.

Mezginejad

Mohammadi

M.H.

Khadem

and Farsani

M.A.

, Evaluation of LKB1 and Serine-Glycine Metabolism Pathway Genes (SHMT1 and GLDC) Expression in AML, Indian J Hematol Blood Transfus 37(2) (2021), 249–255.

65.

Pai

Y.J.

et al., Glycine decarboxylase deficiency causes neural tube defects and features of non-ketotic hyperglycinemia in mice, Nat Commun 6 (2015), 6388.

66.

Mukha

et al., Glycine decarboxylase maintains mitochondrial protein lipoylation to support tumor growth, Cell Metab 34(5) (2022), 775–782e779.

67.

Zhang

W.C.

et al., Glycine decarboxylase activity drives non-small cell lung cancer tumor-initiating cells and tumorigenesis, Cell 148(1-2) (2012), 259–272.

68.

Yuan

et al., Identification of a new GLDC gene alternative splicing variant and its protumorigenic roles in lung cancer, Future Oncol 15(36) (2019), 4127–4139.

69.

Alptekin

et al., Glycine decarboxylase is a transcriptional target of MYCN required for neuroblastoma cell proliferation and tumorigenicity, Oncogene 38(50) (2019), 7504–7520.

70.

Robey

R.W.

Pluchino

K.M.

Hall

M.D.

Fojo

A.T.

Bates

S.E.

and Gottesman

M.M.

, Revisiting the role of ABC transporters in multidrug-resistant cancer, Nat Rev Cancer 18(7) (2018), 452–464.

71.

Locher

K.P.

, Mechanistic diversity in ATP-binding cassette (ABC) transporters, Nat Struct Mol Biol 23(6) (2016), 487–493.

72.

D.M.

Huynh

Truong

A.M.

Haber

and Norris

M.D.

, ABC transporters and neuroblastoma, Adv Cancer Res 125 (2015), 139–170.

73.

Low

F.G.

Shabir

Brown

J.E.

Bill

R.M.

and Rothnie

A.J.

, Roles of ABCC1 and ABCC4 in Proliferation and Migration of Breast Cancer Cell Lines, Int J Mol Sci 21(20) (2020).

74.

Zhao

Guo

Yue

Zhang

and Wang

, ABCC4 is required for cell proliferation and tumorigenesis in non-small cell lung cancer, Onco Targets Ther 7 (2014), 343–351.

75.

Jung

et al., ABCC4/MRP4 contributes to the aggressiveness of Myc-associated epithelial ovarian cancer, Int J Cancer 147(8) (2020), 2225–2238.

76.

Sahores

et al., Multidrug transporter MRP4/ABCC4 as a key determinant of pancreatic cancer aggressiveness, Sci Rep 10(1) (2020), 14217.

77.

Norris

M.D.

et al., Expression of multidrug transporter MRP4/ABCC4 is a marker of poor prognosis in neuroblastoma and confers resistance to irinotecan in vitro , Mol Cancer Ther 4(4) (2005), 547–553.

78.

Murray

et al., Suppression of the ATP-binding cassette transporter ABCC4 impairs neuroblastoma tumour growth and sensitises to irinotecan in vivo , Eur J Cancer 83 (2017), 132–141.

79.

Henderson

M.J.

et al., ABCC multidrug transporters in childhood neuroblastoma: clinical and biological effects independent of cytotoxic drug efflux, J Natl Cancer Inst 103(16) (2011), 1236–1251.

80.

Saygin

Matei

Majeti

Reizes

and Lathia

J.D.

, Targeting Cancer Stemness in the Clinic: From Hype to Hope, Cell Stem Cell 24(1) (2019), 25–40.

81.

Costoya

J.A.

and Arce

V.M.

, Cancer cells escape the immune system by increasing stemness through epigenetic reprogramming, Cell Mol Immunol 20(1) (2022), 6–7.

Machine learning methods revealed the roles of immune-metabolism related genes in immune infiltration,stemness,and prognosis of neuroblastoma

Abstract

BACKGROUND:

METHODS:

RESULTS:

CONCLUSIONS:

Keywords

1. Introduction

2. Materials and methods

2.1 Data acquisition

2.2 Differentially expressed analysis

2.3 Enrichment analysis

2.4 Screening of prognostic genes

2.5 Identification of prognostic gene feature

2.6 Construction of prognosis model

2.7 Validation of prognostic models

2.8 Immune infiltration and stemness analysis

2.9 Correlation analysis of IMRGs with immune score and stemness index

2.10 Statistical analysis

3. Results

3.1 Study cohort

Table 1 Characteristics of study cohort

3.4 Screening of prognostic genes

3.7 Screening of key genes for immune infiltration and stemness

5. Conclusion

Abbreviations

Data and code availability

Author contributions

Funding

Ethics approval and consent to participate

Consent for publication

Competing interests

Supplementary data

sj-xlsx-1-cbm-10.3233_CBM-230119.xlsx - Supplemental material

sj-xlsx-2-cbm-10.3233_CBM-230119.xlsx - Supplemental material

Footnotes

Acknowledgments

References

Table 1
Characteristics of study cohort