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

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:

Risk Score = ∑ 𝐶𝑜𝑒𝑓 ⁢ ( i ) × 𝐸𝑥𝑝𝑟 ⁢ ( i )

𝐶𝑜𝑒𝑓 ⁢ ( i ) refers to coefficient representing each gene obtained by LASSO regression; 𝐸𝑥𝑝𝑟 ⁢ ( i ) 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.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].

OPEN IN VIEWER

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.

OPEN IN VIEWER

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

OPEN IN VIEWER

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).

OPEN IN VIEWER

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.

OPEN IN VIEWER

OPEN IN VIEWER

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 λ 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).

OPEN IN VIEWER

OPEN IN VIEWER

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

OPEN IN VIEWER

OPEN IN VIEWER

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).