Metabolism and Tumor Microenvironment-Associated Molecular Classification of Neuroblastoma

doi:10.21203/rs.3.rs-2748220/v1

Download PDF

Research Article

Metabolism and Tumor Microenvironment-Associated Molecular Classification of Neuroblastoma

https://doi.org/10.21203/rs.3.rs-2748220/v1

This work is licensed under a CC BY 4.0 License

Version 1

posted

You are reading this latest preprint version

Neuroblastoma (NBL) is clinically and biologically heterogeneous, and novel therapies are desperately needed as the poor prognosis of high-risk NBL cases. Increasingly, studies about metabolic reprogramming and tumor microenvironment (TME) open the way to change cancer risk stratification and treatment. Through machine learning, this study identified two metabolic clusters in NBL, which have distinct clinical features, Hallmark pathways and TME. By Weighted gene co-expression network analysis (WGCNA) and Cytoscape, we discovered that RNA and glycosphingolipid (GSL) metabolism play a crucial role in metabolic subtyping. Subsequently, we constructed and verified a risk signature based on key module genes, which performs a good prediction of NBL prognosis. Two risk groups, divided by the median value, are closely associated with clinical features, and risk scores show an inverse correlation with immune infiltration. In combination with gene set enrichment analysis (GSEA) results, we conclude that ribosome biogenesis regulated by rRNA metabolism might be a target for MYCN-amplification NBLs, and GSL metabolism might contribute to TME formation in NBL. Finally, we tried to predict different risk groups' immunotherapy and chemotherapy sensitivity and screened potential targets that might be useful against NBL. In summary, we used multiple bioinformatics analyses to explore the interaction between metabolic processes and TME in NBL and provide new ideas for developing new therapies.

Neuroblastoma

metabolism

TME

Machine Learning

WGCNA

immunotherapy

chemotherapy

Neuroblastoma (NBL) is children's most frequent extracranial solid malignancy (Matthay et al. 2016). It represents 8% of new diagnoses of childhood tumors, accounting for 15% of childhood cancer deaths (Louis and Shohet 2015). NBL is notable for its clinical heterogeneity, ranging from spontaneous regression to widespread metastases (Tolbert and Matthay 2018). The current risk classification systems classify patients into different risk groups based on several prognostic risk features. According to the risk systems, tailored treatments have improved the survival of low- and intermediate-risk patients. However, the long-term survival rate for high-risk patients is less than 50%(Pinto et al. 2015), and the survivors are at risk of severe illnesses due to the intense treatment regimen (Butler et al. 2021). Unfortunately, high-risk NBL accounts for almost 50% of all cases owing to metastatic disease and MYCN amplification(Campbell et al. 2017). Therefore, it is critical to explore novel biomarkers to improve prediction, optimize prognosis and limit exposure to treatment-related toxicities.

Metabolic reprogramming is well-established as a hallmark of cancer (Faubert et al. 2020). Besides, the tumor microenvironment (TME) has been shown to play a critical role in the evolution of tumors (Blavier et al. 2020). In recent years, multiple therapies have been developed to direct the various components of the TME. Meanwhile, there is growing appreciation that tumor progression can be modulated by the metabolism of non-cancer-cell components in TME (Elia and Haigis 2021). However, the interaction of metabolism and TME in tumor progression remains to be fully elucidated.

In this study, we classified NBL patients into two metabolic subtypes. Each was associated with distinct clinical characteristics. Then co-express networks among metabolic processes, Hallmark pathways, and immune signatures were constructed and visualized using WGCNA and Cytoscape. These networks have made significant breakthroughs in our knowledge regarding the TME metabolic process of NBL. Additionally, we developed an MYCN-related metabolic prognostic model and verified the efficiency with an external cohort. Finally, we analyzed the chemotherapy and immunotherapy responses between the two risk groups. These discoveries can effectively predict the prognosis of NBL and provide promising therapies that target high-risk NBL.

2.1 Cohort dataset Download and processing

Neuroblastoma RNA-seq data of GDC TARGET-NBL and GSE62564, together with corresponding clinical information, were obtained from the UCSC Genome Browser database (https://genome.ucsc.edu) and Gene Expression Omnibus (GEO) (https://www.ncbi.nlm.nih.gov/geo/). All sequencing data were converted into transcripts per million (TPM). Then averaged expression values of the same gene were utilized for subsequent studies.

Subsequently, we filtered all samples with the following criteria:
(1) samples with overall survival (OS) time<30 days or incomplete clinical data were excluded.
(2) patients at stages I and II were excluded because of their ideal outcomes.

Finally, 150 samples from the TARGET-NBL dataset were applied as a discovery cohort, and 295 samples from GSE62564 were used as a validation cohort.

2.2 Gene Set Variation Analysis (GSVA)

GO biological processes (GO.BP) gene sets (c5.go.bp.v7.5.1.symbols.gmt) and Hallmark gene sets (h.all.v7.5.1.symbols.gmt) were downloaded from MSigDB v 7.5.1(http://software.broadinstitute.org/gsea/msigdb). A total of 463 metabolic-related gene sets were retrieved from GO.BP gene sets by searching "METABOLIC_PROCESS." GSVA was used to transform the gene expression matrix into a pathway expression matrix by the R package "GSVA."

2.3 Identification and validation of metabolic subtypes

We conducted a Univariable Cox regression analysis using the "survival" package in R software (P<0.05). Accordingly, we identified prognosis-related metabolic gene sets in the TARGET-NBL cohort. Subsequently, the "ConsensusClusterPlus" package in R was used for unsupervised consensus clustering to identify the metabolic classification of the NBL. We employed a defined subsampling parameter of 80% with 100 iterations and set the number of possible clusters (k) to range from 2 to 10 to avoid excessive numbers of groups. The consensus heatmap and cumulative distribution function (CDF) were applied to assess the optimal K. Then, the subtype assignments were validated using T-distributed stochastic neighbor embedding (t-SNE). Valuable clinical data and survival analysis were compared between the various subtypes. Finally, the above approach was also applied to GSE62564 by using the same candidate gene sets.

2.4 Estimation of Tumor Microenvironment (TME) and Tumor Immunity

In the ARGET-NBL cohort, we used the ESTIMATE method to calculate samples' stromal and immune scores. Then the ESTIMATE and consensus measurement of purity estimations (CPE) algorithms were used to calculate the tumor purity (Yoshihara et al. 2013). The above operations were completed through the R package "estimate." In addition, the relative infiltration levels of 23 immune cell types and 13 immunological functions were evaluated by the ssGSEA algorithm in the R package "GSVA." Meanwhile, other tumor immune signatures, including the expression of human leukocyte antigen (HLA) family genes and immune checkpoints, were assessed in NBL samples.

In addition, we use TIMER (https://cistrome.shinyapps.io/timer/) to analyze the immune cell infiltration of samples calculated by different methods (Li et al. 2017).

2.5 Weighted Gene Co-expression Network Analysis (WGCNA)

The R package WGCNA (v1.61) was applied to construct the co-expression network from the TARGET-NBL TPM data (Langfelder and Horvath 2008). In our study, genes in the top quartile of variance and mean expression>1 and the different Hallmark and immune gene sets between two subtypes obtained previously were employed to perform WGCNA. First, the soft-threshold power (β) was calculated based on the scale-free topology fit index (>0.9) to construct a scale-free network. Next, the data matrix was converted into a topological overlap matrix (TOM). All genes and gene sets were classified into distinct modules using the TOM-based dissimilarity measure. Subsequently, we calculated module eigengenes (MEs) and clustered modules according to the MEs' dissimilarity. The following settings were made for the hierarchical clustering:minModuleSize=100; MEDissThres = 0.35. Clinical characteristics, including age, event, fustat, MYCN state, and cluster, were considered traits and corresponded with MEs.P-values and Pearson's correlation test were applied to identify key modules highly relevant to the clinical trait. Co-expression networks(weights >0.1) related to metabolic progress, hallmark pathways, and immune features were selected from key modules. Ultimately, we visualized co-expression relationships between gene expression, signaling pathways, immune signatures, and metabolic processes by Cytoscape. Finally, the key module genes were extracted to further model construction.

2.6 Construction and Verification of the gene signature

First, Univariable Cox regression analysis was performed to identify genes associated with prognosis in the TARGET-NBL. Subsequently, we use the least absolute shrinkage and selection operator (LASSO) Cox regression analysis to minimize overfitting. Then screened genes were employed to construct the prognostic risk signature. The risk score for each patient was calculated by the expression of genes multiplied by the relative coefficient (β) (Risk score = ∑ EXP (mRNA) * β). All patients were classified into high-risk and low-risk subgroups using the median value. In addition, the same formula and cut-off value were applied to the GSE62564 cohorts. Finally, the time-dependent ROC(Receiver Operating Characteristic) curves, Kaplan–Meier survival curves, and log-rank tests were used to assess the performance of the signature in both TARGET-NBL and GSE62564.

2.7 Gene Set Enrichment Analysis(GSEA)

We downloaded The GSEA software from Molecular Signatures Database (http://www.gsea-msigdb.org/gsea/downloads.jsp). Using KEGG gene sets (c2.cp.kegg.v2022.1.Hs.symbols.gmt), GSEA analysis was conducted to investigate differences in gene enrichment between high-risk and low-risk groups.

2.8 Prediction of Immunotherapeutic and Chemotherapy Sensitivity

To further explore new treatments for NBL, Tumor Immune Dysfunction and Exclusion (TIDE) was employed to predict the response of different risk groups to immune checkpoint inhibitors(ICIs)(Jiang et al. 2018b). Each sample got a TIDE score through the TIDE algorithm. A Higher tumor TIDE prediction score is associated with worse ICB response, and patients with a TIDE value <0 were defined as responding to ICIs.

The Connectivity map (https://clue.io/), which allows users to predict compounds that can activate or inhibit based on a gene expression signature, was used to determine which small molecule drugs might be helpful against NBLs(Subramanian et al. 2017). First, DEGs between the high-risk and low-risk groups were obtained with FDR < 0.05 and | log2fold change (FC) | > 1. Subsequently, favorable and unfavorable genes were screened via univariate-cox regression analysis. Furthermore, Two sets of genes above (favorable and unfavorable) were used to acquire candidate drug signatures. Mechanism of actions (MoA) was analyzed based on compounds scoring below -98 in unfavorable gene signature (UGS). Finally, we intersected compounds with a score of <-95 in UGS and compounds of >95 in the favorable gene signature (FGS). These compounds used in clinical therapy are considered candidates for NBL treatment.

To verify the sensitivity of candidate agents, we employed the R package "pRRophetic" to calculate the IC50 of samples in different risk groups by ridge regression(Geeleher et al. 2014). The Spearman's correlation coefficients between drug sensitivity and risk score of each sample were calculated, and P<0.001 was considered statistically significant.

2.9 Statistical Analysis

This study performed all statistical analyses and figures using R (version 4.1.2). The differences between the two groups were compared by using the Wilcoxon test. Kaplan–Meier curves were used to analyze the OS differences and compared using the log-rank test. The univariate Cox and LASSO regression analyses were conducted to identify prognostic factors and metabolic risk models. The predictive performance of models was evaluated through ROC. Pearson's correlation tests were analyzed to calculate the coexpression of different factors. All tests were two‐sided, and P < 0.05 was considered statistically significant.

3.1 Identification and Validation of Metabolism-Associated Subtypes in NBL

In the TARGET-NBL cohort,56 metabolic genesets associated with prognosis were obtained by univariate Cox analysis. The "ConsensusClusterPlus" package in R was used to classify metabolic subtypes of NBL. The cluster effect was best when 150 samples were divided into two subtypes(Figure 1A). The t-SNE showed that most Cluster A and Cluster B samples could be well distinguished, indicating that this subgroup clustering was a reasonable choice (Figure 1B). Kaplan-Meier survival curves were performed to analyze the prognostic characteristics of metabolic subtypes. As shown in (Figure 1C-D), Cluster A presented a more favorable OS and EFS than Cluster B. The same metabolic genesets and cluster methods were employed for molecular subtyping of GSE62564, and similar conclusions were obtained (Figure 1E-H). The heatmap revealed the expression of prognosis-related metabolic genesets in different clusters and the correlation between molecular subtypes and clinical features, such as Age, INSS, and MYCN. As shown in (Figure 1I-J), patients with MYCN amplification, age >1.5 years, and stage IV are primarily concentrated in Cluster B. Meanwhile, metabolic profiles in the different subtypes were consistent between TARGET and GSE62564.

3.2 Hallmark pathway and immune-related characteristics of subtypes

In TARGET-NBL, we performed GSVA to identify Hallmark pathways' differences between two metabolic-related subtypes. As shown in (Figure 2A), the higher correlations of Mitotic Spindle, DNA Repair, G2M Checkpoint, E2F Targets, MYC Targets-V1, MYC Targets-V2, and Spermatogenesis were found in ClusterB compared to ClusterA. On the other hand, other pathways were upregulated in Cluster A. Tumor microenvironment (TME) in NBL samples was estimated via the ESTIMATE algorithm. Our results showed that the immune score, stromal score, and ESTIMATE score displayed higher in Cluster A, and in contrast, tumor purity was higher in Cluster B (Figure 2B). Subsequently, we analyzed the immune features of two subgroups. According to (Figure 2C), most of the immune cells had higher scores in Cluster A, and no significant difference in Activated CD4 T cell and Type 2 T helper cell were detected between the two clusters. Similarly, various immune processes of Cluster A are more vigorous than Cluster B (Figure 2D). Furthermore, it was also shown that most immune checkpoints and HLAs were highly expressed in Cluster A compared to Cluster B (Figure 2E-F). To verify the above conclusions, we got similar results from GSE62564 through the same methods(S1A-F).

3.3 Construction of Networks Among Metabolic processes, Hallmark Pathways, and Immune Signatures

To explore the connection among metabolic processes, Hallmark pathways, and immune signatures, we performed WGCNA on TARGET to construct the co-expression networks. After setting soft-threshold power(β) = 7 (Figure 3A), a scale-free network was constructed based on the scale-free topology R²= 0.92 (Figure S2). Then seven co-expression modules were eventually identified by cutting MEs' clustering height = 0.35 (Figure 3B-C). According to the heatmap of module-trait relationships, the yellow and black module demonstrates the highest correlations with metabolic clusters and MYCN status (Figure 3D). Correlation analysis was performed between gene significance (GS) for cluster and module membership (MM) for genes in each module to test the relationships between the module and metabolic cluster. The results showed that the yellow module (correlation coefficient = 0.63, P value = 1e-48) and the black module(correlation coefficient = 0.68, P value = 2e-172) exhibited the most significant positive correlation with the metabolic cluster (Figure S3-4). Finally, we applied Cytoscape to visualize the networks in the yellow and black modules. It was shown that the genes in the yellow module are mainly involved in the MYCN signaling pathway and the metabolic processes of RNA, pyrimidine, and folic acid (Figure 3E). On the other hand, the genes in the black module are linked to various immune-mediated signaling pathways and metabolic processes (Figure 3F).

3.4 Development of Metabolic Risk Score for NBL

We applied univariate-cox regression analysis to assess the prognostic implications of the key module genes (KMGs). As a result, 18 prognostic-related genes were found in the TARGET dataset (Figure S5). Subsequently, LASSO-Cox regression analysis was utilized to screen the candidate prognostic KMGs (Figure 4A). As a result, four candidate metabolic-related genes were identified for model construction. The risk score formula was as follows: Risk score = 0.0537897200512298 mRNA expression of POLR3D + 0.173843171484948 mRNA expression of PAM16 + 0.138278386088605 mRNA expression of RPL18A + (-0.636625110125905) mRNA expression of ITPRID2. All patients were stratified into high-risk and low-risk subgroups according to median risk scores. Both overall survival (OS) and Event Free survival (EFS) analyses showed that the low-risk group displayed a better prognosis than the high-risk group (Figure 4B). Then, ROC curves were conducted to evaluate the metabolic risk score accuracy and sensitivity. Risk score showed a better predictive performance than Age, INSS, and MYCN. Meanwhile, the AUC values of 1, 3, and 5-year OS were 0.772, 0.780, and 0.713, respectively (Figure 4C). The distribution of survival status and the expression of candidate genes in the groups were visualized in (Figure 4D). It was found that the high-risk group had a significantly increased number of dead status than the low-risk group. In addition, POLR3D, PAM16, and RPL18A were significantly upregulated in high-risk samples, and ITPRID2 was higher expressed in low-risk samples. Finally, the same results were found in GSE62564 (Figure 4E-G).

3.5 Correlation of Risk Model with Clinical and Immune Features

We analyzed the clinical characteristics of different risk groups. The alluvial diagram clearly showed the relationships among the risk groups, metabolic clusters, and patients' overall survival status. Most of the patients in the low-risk group belong to Cluster A, and they were primarily associated with a better prognosis (Figure 5A). Most patients under 1.5 years, III or IVS stage, were enriched in the low-risk group. In contrast, MYCN-amplified patients were more likely to be in the high-risk group(Figure 5B). As shown in (Figure 5C), our model had a good prediction performance in both patients with normal or amplified MYCN. Through GSEA, we found that the high-risk group was associated with cell cycle pathways such as base excision repair, DNA replication, pyrimidine metabolism, ribosome, and RNA polymerase. While neural function pathways(axon guidance, dorsal-ventral axis formation, and long-term potentiation) and immune pathways(FcγR mediated phagocytosis, natural killer cell-mediated cytotoxicity) in the low-risk group were more potent than high-risk group (Figure 5D). Subsequently, We analyzed the immune features of different risk groups. It was confirmed that more immune cell infiltration appeared in the low-risk group of patients using different calculation methods (Figure 5E). Furthermore, Risk scores were inversely proportional to the degree of infiltration of various immune cells (Figure 5F).

3.6 Prediction of Immunotherapeutic and Chemotherapy Sensitivity

In TARGET-NBL, we evaluated the expression of immune checkpoints between the two prognostic risk subgroups. As shown in (Figure 6A), most of the checkpoints had higher expression in the low-risk group. According to the TIDE algorithm, almost all patients in both high-risk and low-risk groups exhibited a positive response to ICIs. Moreover, the TIDE score was lower in the high-risk group (Figure 6B). Overall, ICIs are effective in most NBL patients, and the high-risk group is more sensitive to immunotherapy. Based on the connectivity map (CMap), mechanism of action (MoA) analysis was employed to identify potential compounds capable of targeting high-risk groups. The mechanisms of action, which contained >3 compounds, were screened, and we found 12 targets shared by 55 compounds.14 compounds shared the MoA of HDAC inhibitor. Aurora kinase, Topoisomerase, CDK, FLT3, MEK, and VEGFR contained 6 compounds as the target (Figure 6C). Subsequently, we analyzed the CMap score of drugs commonly used in NBL and pediatric malignant tumors (Figure S6-7). Furthermore, we investigated drugs used in clinical therapy and intersected drugs with a score of <-95 in UGS and drugs of >95 in the favorable gene signature (FGS) (Figure 6D) and finally, employed the R package "pRRophetic" to verify the sensitivity of candidate drugs. The results demonstrated the clinical implications of these drugs and were consistent with the conclusions of CMap (Figure 6E).

Neuroblastoma is children's most common extracranial solid tumor, responsible for approximately 15% of all childhood cancer deaths. Although the current risk stratifications and therapies have improved the prognosis of low-risk and intermediate-risk NBL, the outcome of high-risk NBL is still dismal. Thus, exploring solutions for improving patient survival and reducing treatment‑associated toxicities is essential. In addition, metabolic reprogramming of organisms, particularly non-cancer-cell types in TIME, has been confirmed to play a crucial role in tumor initiation and progression. However, the relationships between metabolism and TME in NBL remain unclear.

In this study, we classified NBLs into two subtypes based on unsupervised consensus clustering analysis and found significant differences in the clinical characteristics of these two subtypes. Almost all patients harboring MYCN amplification were categorized into Cluster B, which indicates that MYCN is a decisive gene regulating the metabolic reprogramming of NBL. MYCN is an oncogenic transcription factor that regulates the transcription of many genes involved in fundamental cellular processes, including metabolism, cell growth, proliferation, differentiation, and migration (Ruiz-Pérez et al. 2017). Of note, MYCN amplification NBLs, accounting for 20% of all cases, represent the most malignant NBL subtype (Smith and Foster 2018). It was previously reported that E2F and MYC are functionally associated (Strieder and Lutz 2003); DNA repair and G2M checkpoint are correlated with MYCN (Hu et al. 2019; Aygun and Altungoz 2019). Furthermore, MYCN is essential in remodeling the immunological microenvironment toward a suppressive phenotype (Raieli et al. 2021). These results are consistent with our results.

Subsequently, we performed WGCNA to construct the co-expression networks among metabolic processes, Hallmark pathways, and immune signatures. Our results revealed that the RNA (particularly rRNA) metabolism- and cell cycle-related metabolic processes, driven by MYCN, play a crucial role in Cluster B. In contrast, metabolic processes related to glycosphingolipid (GSL) and glycosylceramide (GC) participate in the tumor microenvironment (TME) formation in Cluster A through various pathways. Moreover, it is well-documented that RNA metabolisms, transport, degradation, and microRNA processing, are critical steps of gene expression (Rocak and Linder 2004). Ribosomal RNAs (rRNAs) and transfer RNAs (tRNAs) are essential components of ribosome biogenesis and polymerase-regulated gene transcription (Johnson and Johnson 2008). Ribosome biogenesis is essential in the tumorigenic process and is frequently associated with poor prognosis(Pelletier et al. 2018). It is well-established that MYC is the master regulator of ribosome biogenesis and controls rRNA processing (van Riggelen et al. 2010; Destefanis et al. 2020). Thus, we believe that ribosome biogenesis, regulated by MYCN via modulating rRNA metabolism, is an important foundation for metabolic typing. Glycosphingolipids, normally found at the outer surface of cell membranes, are believed to play an important role in receptor function and cell immune recognition (Lalazar et al. 2006). However, the causal connection between the glycosphingolipid and TME in NBL has been elusive. Our studies show a possibility that the TME of NBL can be improved by regulating GSL metabolism.

Then, we identified four genes, including POLR3D, PAM16, RPL18A, and ITPRID2, for constructing a LASSO-cox model to predict clinical outcomes in NBL patients. This model showed good predictive performance for NBL prognostic and was validated by the GSE62564 cohorts. Based on the median value of the risk score, we divided all patients into high-risk and low-risk subgroups and compared their manifestations, signaling pathways, and immune cell infiltration between different subtypes. We found that these two subgroups are closely related to clinical features, including outcome, age, INSS, and MYCN. In addition, risk scores are inversely proportional to the degree of infiltration of various immune systems. Through GSEA, we found that the high-risk group is associated with gene expression pathways such as pyrimidine metabolism, base excision repair, DNA replication, RNA polymerase, and the ribosome. These results support our hypothesis that ribosome biogenesis is an important foundation for metabolic typing. In contrast, neurodevelopmental pathways, such as axon guidance, dorsoventral axis formation, and long-term potentiation, combine with immune regulatory pathways, including FCγR mediated phagocytosis and natural killer (NK) cell-mediated cytotoxicity. Thus, we believe GSLs modulate TME in NBL via pathways related to FCγR mediated phagocytosis and NK cell-mediated cytotoxicity.

POLR3D, also known as BN51, is a co-factor of RNA polymerase III and is involved in rRNA metabolism (Greasley et al. 2000). RPL18A is a 60S ribosomal subunit protein family member involved in rRNA processing steps and ribosome biogenesis (Leh et al. 2000; Konikkat and Woolford 2017). Our findings through WGCNA were consistent with these reports. Thus, our method of constructing metabolic-immune network practicability was verified. PAM16 (Presequence Translocase Associated Motor 16) is also known as MAGMAS(Mitochondria-Associated Granulocyte Macrophage CSF-Signaling Molecule) and is involved in mitochondrial protein transport (Jubinsky et al. 2001). The expression pattern of Magmas has been reported to be associated with chemotherapy resistance in OC cells(Ahmed et al. 2020), but little is known about NBL. ITPRID2, also named sperm-specific antigen 2 (SSFA2), was reported to be related to calcium signaling and the EMT pathway (Okayama et al. 2016), and its role in human tumors remains unclear.

Furthermore, we employed TIDE to predict the response of different risk groups to ICIs and used CMap to screen drugs that might be useful against NBLs. We found that almost all patients might respond positively to ICIs. Therefore, drugs targeting HDAC, Aurora, CDK, FLT3, MEK, and Topoisomerase might be treatment options for NBL. In this study, we confirmed the effectiveness and limitations of the present drugs applied for NBL and provided candidate drugs that might be useful for NBL. These candidate drugs are supported by clinical trials or experimental evidence (Minturn et al. 2011; McFarland et al. 2018; Jiang et al. 2018a; DuBois et al. 2018; He et al. 2021; Foster et al. 2021; Dedoni et al. 2021).

However, there are certain deficiencies in our study. First, we performed research based on bioinformatics analyses, and further validations in vivo and in vitro are required. Besides, more samples and clinical trials are needed to confirm our results. Finally, we will verify the possible mechanism of the metabolic-TME network in the future.

Funding

This work was funded by the Shanghai Shen Kang Hospital Development Ccenter (Grant numbers: SHDC22020210) and Shanghai Children's Hospital(2020YGZM05). Sponsoring foundations exert no effect on the writing of this manuscript.

Competing Interests

The authors have no relevant financial or non-financial interests to disclose.

Author Contributions

Xin Du designed the study, formed bioinformatics analysis and wrote the manuscript. JB S,SY J, JW Y and XL L provided clinical guidance and made the literature search. Bing Zou prepared the figures.

Data Availability

Ethics approval

Expression and clinical data are from the database, and there are no ethical issues.

Consent to participate

The authors affirm data in this manuscript are from databases. This research is not involving human subjects.

Consent to publish

The authors affirm data in this manuscript are from databases and do not contain any person’s data.

Ahmed N, Kadife E, Raza A, et al (2020) Ovarian Cancer, Cancer Stem Cells and Current Treatment Strategies: A Potential Role of Magmas in the Current Treatment Methods. Cells 9:E719. https://doi.org/10.3390/cells9030719
Aygun N, Altungoz O (2019) MYCN is amplified during S phase, and c–myb is involved in controlling MYCN expression and amplification in MYCN–amplified neuroblastoma cell lines. Mol Med Rep 19:345–361. https://doi.org/10.3892/mmr.2018.9686
Blavier L, Yang R-M, DeClerck YA (2020) The Tumor Microenvironment in Neuroblastoma: New Players, New Mechanisms of Interaction and New Perspectives. Cancers (Basel) 12:E2912. https://doi.org/10.3390/cancers12102912
Butler E, Ludwig K, Pacenta HL, et al (2021) Recent progress in the treatment of cancer in children. CA Cancer J Clin 71:315–332. https://doi.org/10.3322/caac.21665
Campbell K, Gastier-Foster JM, Mann M, et al (2017) Association of MYCN copy number with clinical features, tumor biology, and outcomes in neuroblastoma: A report from the Children’s Oncology Group. Cancer 123:4224–4235. https://doi.org/10.1002/cncr.30873
Dedoni S, Marras L, Olianas MC, et al (2021) The Neurotrophin Receptor TrkC as a Novel Molecular Target of the Antineuroblastoma Action of Valproic Acid. Int J Mol Sci 22:7790. https://doi.org/10.3390/ijms22157790
Destefanis F, Manara V, Bellosta P (2020) Myc as a Regulator of Ribosome Biogenesis and Cell Competition: A Link to Cancer. Int J Mol Sci 21:E4037. https://doi.org/10.3390/ijms21114037
DuBois SG, Mosse YP, Fox E, et al (2018) Phase II Trial of Alisertib in Combination with Irinotecan and Temozolomide for Patients with Relapsed or Refractory Neuroblastoma. Clin Cancer Res 24:6142–6149. https://doi.org/10.1158/1078-0432.CCR-18-1381
Elia I, Haigis MC (2021) Metabolites and the tumour microenvironment: from cellular mechanisms to systemic metabolism. Nat Metab 3:21–32. https://doi.org/10.1038/s42255-020-00317-z
Faubert B, Solmonson A, DeBerardinis RJ (2020) Metabolic reprogramming and cancer progression. Science 368:eaaw5473. https://doi.org/10.1126/science.aaw5473
Foster JH, Voss SD, Hall DC, et al (2021) Activity of Crizotinib in Patients with ALK-Aberrant Relapsed/Refractory Neuroblastoma: A Children’s Oncology Group Study (ADVL0912). Clin Cancer Res 27:3543–3548. https://doi.org/10.1158/1078-0432.CCR-20-4224
Geeleher P, Cox N, Huang RS (2014) pRRophetic: an R package for prediction of clinical chemotherapeutic response from tumor gene expression levels. PLoS One 9:e107468. https://doi.org/10.1371/journal.pone.0107468
Greasley PJ, Bonnard C, Amati B (2000) Myc induces the nucleolin and BN51 genes: possible implications in ribosome biogenesis. Nucleic Acids Res 28:446–453. https://doi.org/10.1093/nar/28.2.446
He D, Fu S, Zhou A, et al (2021) Camptothecin Regulates Microglia Polarization and Exerts Neuroprotective Effects via Activating AKT/Nrf2/HO-1 and Inhibiting NF-κB Pathways In Vivo and In Vitro. Front Immunol 12:619761. https://doi.org/10.3389/fimmu.2021.619761
Hu X, Zheng W, Zhu Q, et al (2019) Increase in DNA Damage by MYCN Knockdown Through Regulating Nucleosome Organization and Chromatin State in Neuroblastoma. Front Genet 10:684. https://doi.org/10.3389/fgene.2019.00684
Jiang J, Song X, Yang J, et al (2018a) Triptolide Inhibits Proliferation and Migration of Human Neuroblastoma SH-SY5Y Cells by Upregulating MicroRNA-181a. Oncol Res 26:1235–1243. https://doi.org/10.3727/096504018X15179661552702
Jiang P, Gu S, Pan D, et al (2018b) Signatures of T cell dysfunction and exclusion predict cancer immunotherapy response. Nat Med 24:1550–1558. https://doi.org/10.1038/s41591-018-0136-1
Johnson DL, Johnson SAS (2008) Cell biology. RNA metabolism and oncogenesis. Science 320:461–462. https://doi.org/10.1126/science.1158680
Jubinsky PT, Messer A, Bender J, et al (2001) Identification and characterization of Magmas, a novel mitochondria-associated protein involved in granulocyte-macrophage colony-stimulating factor signal transduction. Exp Hematol 29:1392–1402. https://doi.org/10.1016/s0301-472x(01)00749-4
Konikkat S, Woolford JL (2017) Principles of 60S ribosomal subunit assembly emerging from recent studies in yeast. Biochem J 474:195–214. https://doi.org/10.1042/BCJ20160516
Lalazar G, Preston S, Zigmond E, et al (2006) Glycolipids as immune modulatory tools. Mini Rev Med Chem 6:1249–1253. https://doi.org/10.2174/138955706778742722
Langfelder P, Horvath S (2008) WGCNA: an R package for weighted correlation network analysis. BMC Bioinformatics 9:559. https://doi.org/10.1186/1471-2105-9-559
Leh V, Yot P, Keller M (2000) The cauliflower mosaic virus translational transactivator interacts with the 60S ribosomal subunit protein L18 of Arabidopsis thaliana. Virology 266:1–7. https://doi.org/10.1006/viro.1999.0073
Li T, Fan J, Wang B, et al (2017) TIMER: A Web Server for Comprehensive Analysis of Tumor-Infiltrating Immune Cells. Cancer Res 77:e108–e110. https://doi.org/10.1158/0008-5472.CAN-17-0307
Louis CU, Shohet JM (2015) Neuroblastoma: molecular pathogenesis and therapy. Annu Rev Med 66:49–63. https://doi.org/10.1146/annurev-med-011514-023121
Matthay KK, Maris JM, Schleiermacher G, et al (2016) Neuroblastoma. Nat Rev Dis Primers 2:16078. https://doi.org/10.1038/nrdp.2016.78
McFarland AJ, Davey AK, McDermott CM, et al (2018) Differences in statin associated neuroprotection corresponds with either decreased production of IL-1β or TNF-α in an in vitro model of neuroinflammation-induced neurodegeneration. Toxicol Appl Pharmacol 344:56–73. https://doi.org/10.1016/j.taap.2018.03.005
Minturn JE, Evans AE, Villablanca JG, et al (2011) Phase I trial of lestaurtinib for children with refractory neuroblastoma: a new approaches to neuroblastoma therapy consortium study. Cancer Chemother Pharmacol 68:1057–1065. https://doi.org/10.1007/s00280-011-1581-4
Okayama A, Kimura Y, Miyagi Y, et al (2016) Relationship between phosphorylation of sperm-specific antigen and prognosis of lung adenocarcinoma. J Proteomics 139:60–66. https://doi.org/10.1016/j.jprot.2016.03.005
Pelletier J, Thomas G, Volarević S (2018) Ribosome biogenesis in cancer: new players and therapeutic avenues. Nat Rev Cancer 18:51–63. https://doi.org/10.1038/nrc.2017.104
Pinto NR, Applebaum MA, Volchenboum SL, et al (2015) Advances in Risk Classification and Treatment Strategies for Neuroblastoma. J Clin Oncol 33:3008–3017. https://doi.org/10.1200/JCO.2014.59.4648
Raieli S, Di Renzo D, Lampis S, et al (2021) MYCN Drives a Tumor Immunosuppressive Environment Which Impacts Survival in Neuroblastoma. Front Oncol 11:625207. https://doi.org/10.3389/fonc.2021.625207
Rocak S, Linder P (2004) DEAD-box proteins: the driving forces behind RNA metabolism. Nat Rev Mol Cell Biol 5:232–241. https://doi.org/10.1038/nrm1335
Ruiz-Pérez MV, Henley AB, Arsenian-Henriksson M (2017) The MYCN Protein in Health and Disease. Genes (Basel) 8:E113. https://doi.org/10.3390/genes8040113
Smith V, Foster J (2018) High-Risk Neuroblastoma Treatment Review. Children (Basel) 5:E114. https://doi.org/10.3390/children5090114
Strieder V, Lutz W (2003) E2F proteins regulate MYCN expression in neuroblastomas. J Biol Chem 278:2983–2989. https://doi.org/10.1074/jbc.M207596200
Subramanian A, Narayan R, Corsello SM, et al (2017) A Next Generation Connectivity Map: L1000 Platform and the First 1,000,000 Profiles. Cell 171:1437–1452.e17. https://doi.org/10.1016/j.cell.2017.10.049
Tolbert VP, Matthay KK (2018) Neuroblastoma: clinical and biological approach to risk stratification and treatment. Cell Tissue Res 372:195–209. https://doi.org/10.1007/s00441-018-2821-2
van Riggelen J, Yetil A, Felsher DW (2010) MYC as a regulator of ribosome biogenesis and protein synthesis. Nat Rev Cancer 10:301–309. https://doi.org/10.1038/nrc2819
Yoshihara K, Shahmoradgoli M, Martínez E, et al (2013) Inferring tumour purity and stromal and immune cell admixture from expression data. Nat Commun 4:2612. https://doi.org/10.1038/ncomms3612

No competing interests reported.

S1.tif
Figure S1: Relationship between subtypes, Hallmark pathways, and immune-related characteristics in GSE62564.
S2.tif
Figure S2: Distribution of nodes with the degree of connection(k), and correlation of log (k) and log [P(k)].
S3.tif
Figure S3: Scatterplots of GS for the cluster (y-axis) versus MM (x-axis) in the yellow module.
S4.tif
Figure S4: Scatterplots of GS for the cluster (y-axis) versus MM (x-axis) in the black module.
S5.tif
Figure S5: Forest plot of module-related prognostic genes based on univariate Cox regression analysis.
S6.tif
Figure S6: CMap score of drugs commonly used in NBL.
S7.tif
Figure S7: CMap score of drugs commonly used in pediatric malignant tumors.

Download PDF

Version 1

posted

You are reading this latest preprint version

Metabolism and Tumor Microenvironment-Associated Molecular Classification of Neuroblastoma

Status:

Version 1

Abstract

Figures

1 Introduction

2 Materials And Methods

2.1 Cohort dataset Download and processing

2.2 Gene Set Variation Analysis (GSVA)

2.3 Identification and validation of metabolic subtypes

2.4 Estimation of Tumor Microenvironment (TME) and Tumor Immunity

2.5 Weighted Gene Co-expression Network Analysis (WGCNA)

2.6 Construction and Verification of the gene signature

2.7 Gene Set Enrichment Analysis(GSEA)

2.8 Prediction of Immunotherapeutic and Chemotherapy Sensitivity

2.9 Statistical Analysis

3 Results

3.1 Identification and Validation of Metabolism-Associated Subtypes in NBL

3.2 Hallmark pathway and immune-related characteristics of subtypes

3.3 Construction of Networks Among Metabolic processes, Hallmark Pathways, and Immune Signatures

3.4 Development of Metabolic Risk Score for NBL

3.5 Correlation of Risk Model with Clinical and Immune Features

3.6 Prediction of Immunotherapeutic and Chemotherapy Sensitivity

4 Discussion

Declarations

Funding

Competing Interests

Author Contributions

Data Availability

Ethics approval

Consent to participate

Consent to publish

References

Additional Declarations

Supplementary Files

Status:

Version 1