Transcription factor activity-based stratification reveals prognostic subtypes and immune landscape in neuroblastoma
Highlight box
Key findings
• Two transcription factor (TF)-defined neuroblastoma (NB) subtypes were identified; the MYC/E2F-dominant subtype had a 5-year overall survival of 34% vs. 84% in the favorable group.
• A single, principal component analysis-derived TF_score outperformed stage and MYCN status for survival prediction (5-year receiver operating characteristic =0.83) and remained significant in multivariable analysis.
• TF_score-high tumors displayed lower immune-cell infiltration, higher Tumor Immune Dysfunction and Exclusion scores, and weaker predicted response to checkpoint blockade.
What is known and what is new?
• Age, stage, and MYCN amplification only partially explain NB heterogeneity; TF dysregulation is critical but seldom quantified clinically.
• This study provides the first system-level TF-activity map in NB, introduces a reproducible TF_score, and validates its prognostic and immunological relevance across bulk and single-cell datasets.
What is the implication, and what should change now?
• TF_score can be calculated from routine transcriptomes and incorporated into nomograms to refine the current International Neuroblastoma Risk Group risk groups.
• Prospective implementation of TF_score may guide treatment escalation or de-intensification and stratify patients for immunotherapy or MYC/E2F-targeted trials.
Introduction
Neuroblastoma (NB) has an annual incidence of ~10.5 cases per million children under 15 years, with marked geographic and ethnic variation-rates are highest in North America and Northern Europe, intermediate in East Asia, and lowest in sub-Saharan Africa (1). Across regions, NB accounts for 8% of pediatric cancers yet nearly 15% of cancer-related deaths, reflecting its aggressive biology (2,3). Recent large-scale sequencing efforts have expanded the genomic landscape beyond MYCN amplification: whole-genome analyses identified recurrent alterations in telomere-maintenance pathways, chromatin remodelers, and ALK/RAS-MAPK signaling (4), while integrative multi-omic profiling showed that a subset of high-risk tumors lacks canonical driver mutations but displays pervasive transcriptional dysregulation (5). These data underscore that static genomic markers only partially capture the mechanisms governing tumor behavior and highlight the need for functional metrics—such as transcription factor (TF) activity—to refine risk assessment and guide therapy.
TFs orchestrate gene-expression networks that regulate proliferation, lineage commitment, immune signaling, and metabolic rewiring. Aberrant activity of MYCN, PHOX2B, TEAD-YAP, and GATA families has been implicated in NB pathogenesis (6-9); however, most studies have focused on individual TFs or messenger RNA (mRNA) abundance, which poorly reflects true regulatory activity. Emerging algorithms such as decoupleR infer TF activity by integrating weighted TF-target networks with transcriptome profiles, yielding sample-specific activity scores (10). System-level TF activity landscapes remain unexplored in NB, particularly in relation to the tumor immune microenvironment.
Conventional markers (age, stage, and MYCN) only partially explain outcome variability in NB. We hypothesized that sample-level TF activity captures regulatory programs linked to tumor aggressiveness and immune contexture, and that a composite TF-based score would augment current risk models. Building on this rationale, we sought to infer TF activities with decoupleR and delineate TF activity-defined NB subtypes by consensus clustering in GSE49710; to derive a continuous TF_score via Boruta-guided principal component analysis (PCA), compare its prognostic performance with stage and MYCN, and integrate it into an externally validated dataset (E-MTAB-8248). By uniting TF-activity profiling with clinical variables, the study aims to provide a reproducible and practicable framework for refined risk stratification and to support personalized therapeutic decision-making in NB. We present this article in accordance with the TRIPOD reporting checklist (available at https://tp.amegroups.com/article/view/10.21037/tp-2025-aw-739/rc).
Methods
Data acquisition and preprocessing
Expression matrices and matched clinical information from GSE49710 and E-MTAB-8248 datasets were downloaded from Gene Expression Omnibus (GEO) and ArrayExpress. Data preprocessing involved normalization and annotation using standard chip annotation files. Genes with multiple probes were represented by their average expression. GSE49710 was used as the training dataset, while E-MTAB-8248 served as external validation.
GSE49710 comprises 498 primary NB tumors prospectively registered by the Children’s Oncology Group (COG) and collaborating North American centers between 2004 and 2012. E-MTAB-8248 contains 223 tumors collected by the German NB2004 trial network and affiliated European institutions (Germany, Austria, Spain, and the Netherlands) during 2001–2014. All patients were of European or North American ancestry; treatment protocols followed COG or German Society for Pediatric Oncology and Hematology (GPOH) guidelines, respectively.
Single-cell RNA sequencing (RNA-seq) data (Jansky dataset) were acquired through a dedicated Shiny app, following previously described preprocessing methods for filtering, normalization, and annotation.
The study was conducted in accordance with the Declaration of Helsinki and its subsequent amendments. All datasets used in this study were obtained from publicly available repositories and contained only de-identified data; therefore, separate institutional review board approval and informed consent were not required.
TF activity calculation
We retrieved a comprehensive transcription-factor (TF)-target regulon from OmniPathR. TF activities for every sample in both the bulk (GSE49710) and single-cell (Jansky et al.) datasets were quantified with the decoupleR package using its built-in univariate linear model (ULM) algorithm. Briefly, ULM regresses gene-expression values against signed, weighted TF-target interactions and produces a normalized activity z-score after 1,000 gene-label permutations. To facilitate independent interrogation, we release a complete sample-wise TF-activity matrix (498 tumors) in table available at https://cdn.amegroups.cn/static/public/tp-2025-aw-739-1.xlsx. Throughout, we report TF activity inferred from regulon-target relationships; this is distinct from mRNA expression and better reflects functional regulation.
Unsupervised clustering and TF_score construction
We first performed univariate Cox regression on all TF activities in GSE49710; TFs with P<1×10−4 were retained. Activity scores were median-centered and subjected to consensus clustering with ConsensusClusterPlus (Pearson distance, 1,000 resamplings, pItem =0.8, pFeature =1) (11). The cumulative distribution function curve and ΔK statistic indicated that a two-cluster solution optimally captured TF-driven heterogeneity.
To quantify the transcriptional phenotype of individual tumors, we followed the strategy described by Zhang et al. (12). TFs positively correlated with the high-risk cluster were denoted set A, whereas negatively correlated TFs constituted set B. Boruta (500 iterations, default parameters) removed redundant features within each set. PCA was then applied separately to the Boruta-confirmed TFs of sets A and B, and the first principal component (PC1) of each set was extracted. TF_score was defined as PC1positive − PC1negative. Patients were categorized into TF_score-high and TF_score-low groups using the cohort median as a cut-off, enabling downstream survival and immune analyses.
Immune infiltration and response assessment
The Estimation of STromal and Immune cells in MAlignant Tumors using Expression data (ESTIMATE) algorithm was applied to evaluate stromal and immune cell content (13). Specific immune cell subtypes were characterized using Cell-type Identification By Estimating Relative Subsets of RNA Transcripts (CIBERSORT) and Microenvironment Cell Populations-counter (MCP-counter) (14,15). Potential immunotherapy responses were assessed by Immunophenoscore (IPS) and Tumor Immune Dysfunction and Exclusion (TIDE) algorithms (16,17).
Nomogram development and validation
Significant prognostic factors identified by univariate Cox regression (P<0.05) were integrated into a prognostic nomogram. Receiver operating characteristic (ROC) curves, calibration plots, and decision curve analyses (DCAs) were employed for validation.
Statistical analysis
Statistical analyses utilized R software (version 4.2.2). Differences between groups were assessed using the Wilcoxon, Chi-squared, or Fisher’s exact tests. Kaplan-Meier survival analyses, log-rank tests, and Cox proportional hazards models determined prognostic significance, with P values <0.05 considered significant. For every analysis entailing ≥5 concurrent comparisons, raw P values were adjusted using the Benjamini-Hochberg false discovery rate (FDR) method; results were deemed significant at FDR <0.05.
Results
Identification of NB subtypes based on TF activity
We first calculated the TF activities in the GSE49710 dataset using the decoupleR package and performed unsupervised clustering analysis. The full activity values for all TFs across all samples are provided in table available at https://cdn.amegroups.cn/static/public/tp-2025-aw-739-1.xlsx. Consensus clustering analysis revealed two robust and stable clusters (cluster 1 and cluster 2; Figure 1A,1B). Visualization by t-distributed stochastic neighbor embedding (t-SNE) and PCA demonstrated clear separation between the two clusters (Figure 1C,1D). Kaplan-Meier analysis showed significantly worse overall survival (OS) in cluster 2 [hazard ratio (HR) =3.59; P<0.001; Figure 1E]. Cluster 2 was enriched in MYCN amplification and patients aged ≥18 months, suggesting an aggressive phenotype (Figure 1F). Furthermore, hallmark pathway analysis showed distinct functional profiles between the two clusters, with cluster 2 exhibiting higher mRNA stemness index (mRNAsi) scores, indicating greater tumor cell stemness (Figure 1G,1H).
Clinically, the stratification of NB patients into TF activity-based subtypes offers improved risk identification for aggressive disease, particularly those who may not exhibit MYCN amplification but harbor transcriptional signatures indicating stemness and poor prognosis. This has potential implications for identifying candidates for intensified or novel therapeutic regimens.
Quantification and validation of TF activity-based clustering
To further quantify the clustering results, we used Boruta feature selection for dimension reduction and subsequently calculated a TF_score via PCA (Figure 2A). The distribution of TF_score across age, stage, and MYCN status is shown in Figure 2B. This TF_score effectively discriminated between the two initial clusters, with Cluster 2 having significantly higher TF_scores (P<2.2e−16; Figure 2C). Using the median TF_score as a cutoff, patients were divided into TF_score-high and TF_score-low groups, with the TF_score-high group demonstrating significantly worse OS (HR =4.82; P<0.001; Figure 2D). Time-dependent ROC analysis showed that TF_score consistently outperformed MYCN amplification status and tumor stage as a prognostic indicator, with superior predictive accuracy across various time points (Figure 2E). DCA further supported the superior prognostic decision-making capability of TF_score relative to MYCN and stage (Figure 2F). These findings were validated in subgroup analyses by MYCN status, tumor stage, and age (Figure S1), as well as external validation using the E-MTAB-8248 dataset, where TF_score-high patients consistently showed poorer outcomes (Figure S2).
These findings provide a clinically actionable prognostic biomarker that integrates underlying regulatory mechanisms, offering oncologists a tool to more precisely stratify patients and predict treatment outcomes.
Immune landscape analysis based on TF_score
To investigate differences in immune infiltration and immune response between TF_score groups, we applied the ESTIMATE, CIBERSORT, and MCP-counter algorithms. The TF_score-high group demonstrated lower immune infiltration across multiple cell types (Figure 3A). Additionally, analysis of IPS and TIDE scores indicated that TF_score-high patients exhibited higher immune suppression and less favorable immune responsiveness (higher TIDE, lower IPS; Figure 3B,3C). Consistently, TIDE-based immune response prediction showed significantly lower immune response rates in the TF_score-high group compared to the TF_score-low group (Figure 3D,3E).
The immune analyses have potential translational significance. Patients with high TF_scores may be less likely to benefit from immune checkpoint blockade due to a more immunosuppressive tumor microenvironment. Conversely, low TF_score patients may be better candidates for immunotherapy, which may influence clinical decision-making in future trials.
Construction and validation of a prognostic nomogram
Integrating TF_score with clinical factors, including MYCN amplification status and tumor stage, we constructed a prognostic nomogram to predict patient survival outcomes (Figure 4A). ROC and calibration curve analyses confirmed the predictive reliability of the nomogram (Figure 4B,4C). Stratification using median nomogram scores showed robust survival differentiation, with the high-risk group having markedly worse outcomes (HR =13.86, P<0.001; Figure 4D). Further analysis revealed key TFs distinguishing high-risk from low-risk groups, notably elevated activity of MYC and E2F family TFs and reduced activity of FOXO3 and IRF1 in the high-risk group (Figure 4E). Detailed analyses of TF-target interactions further supported the regulatory roles of these TFs (Figure 4F, Figure S3).
This nomogram could potentially be incorporated into clinical workflows to guide risk-adapted therapeutic strategies and follow-up intensity, thus translating omics data into practical clinical applications.
Single-cell data validation
Single-cell RNA-seq data (Jansky dataset) confirmed the differential activity patterns observed in bulk RNA-seq. Uniform Manifold Approximation and Projection (UMAP) visualization showed clear distinctions in TF activities between high- and low-risk tumor cells, with high-risk cells exhibiting increased MYC activity and decreased FOXO3 activity, consistent with bulk transcriptomic findings (Figure 5A-5D). These results underscore the robustness and biological relevance of TF activity-based classification in NB.
The consistency between bulk and single-cell findings enhances the credibility of the TF_score model and suggests its applicability across cellular hierarchies, further supporting its potential integration into translational and clinical research.
Discussion
This study provides a comprehensive, activity-level view of transcriptional regulation in NB. Using decoupleR, we identified two robust TF-defined subtypes and derived a patient-level TF_score that quantitatively captures this phenotype. TF_score consistently outperformed stage and MYCN for survival prediction, remained significant in multivariable analyses and across clinical subgroups, and generalized to an external cohort. Single-cell data corroborated enrichment of MYC/E2F programs with concomitant suppression of FOXO3/IRF1 in high-risk tumor cells. Integration of TF_score with routine variables into a nomogram yielded accurate 1-, 3-, and 5-year survival estimates with strong calibration.
Mechanistically, the aggressive cluster and TF_score-high tumors were dominated by heightened MYC and E2F activity, partners long recognized for accelerating S-phase entry, replication stress, and metabolic re-wiring. Concordantly, these tumors displayed the highest mRNAsi stemness indices, suggesting that MYC/E2F transcriptional programs sustain an undifferentiated, self-renewing cell state. In contrast, FOXO3 and IRF1—gatekeepers of oxidative—stress responses, type I interferon signaling, and major histocompatibility complex (MHC) class I antigen presentation, were suppressed. The reciprocal regulation of these TF sets provides a plausible explanation for the coupling of proliferation, stemness, and immune evasion observed in high-risk NB. Notably, single-cell data confirmed that MYC-high/FOXO3-low states are enriched specifically in high-risk tumor cells rather than the microenvironment, underscoring their tumor-intrinsic origin.
Immune deconvolution revealed that TF_score-high tumors are immunologically “cold”, harboring fewer cytotoxic T cells and antigen-presenting macrophages, lower IPSs, and higher TIDE dysfunction/exclusion scores. From a translational perspective, our data suggest that patients with low TF_scores may derive durable benefit from checkpoint blockade, whereas those with high TF_scores are unlikely responders unless combination regimens convert the microenvironment to an inflamed state. Pre-clinical studies combining MYC inhibition with STING or CD40 agonists have shown promise in MYC-driven malignancies (18,19), offering a rational avenue for TF_score-guided immunotherapy trials in NB.
Although the present study did not explicitly re-classify patients according to the International Neuroblastoma Risk Group (INRG) criteria, our data suggest that TF_score captures biological risk components not encompassed by age, stage, or MYCN status. Prospective cohorts that record full INRG variables will be needed to determine whether TF_score can re-identify intermediate-risk patients with adverse biology or high-risk patients with favorable biology, thereby refining treatment decisions. TF_score can be computed from standard transcriptomes and incorporated into existing workflows to refine risk stratification, support discussions at multidisciplinary boards, and inform trial allocation-identifying patients more likely to benefit from treatment escalation or, conversely, candidates for de-intensification. Immune profiling linked to TF_score also provides a pragmatic readout to anticipate response to checkpoint blockade.
Hyperdiploid/near-triploid NBs, characterized mainly by whole-chromosome (numeric) gains, generally have favorable outcomes, whereas near-diploid tumors enriched for segmental chromosomal alterations exhibit inferior survival. TF activity, however, reflects network-level regulation and can be increased by copy-number-neutral mechanisms (e.g., super-enhancers, kinase-mediated TF activation, epigenetic remodeling). Our TF_score captures these program-level changes and therefore complements, rather than contradicts, established ploidy-based risk associations.
The TF landscape also nominates therapeutically tractable vulnerabilities. Direct MYC targeting remains challenging, yet indirect blockade via BET bromodomain, CDK7/9, or AURKA inhibitors can down-regulate MYC/E2F transcriptional output and has entered early-phase pediatric trials (20-22). Our data predict that such agents would preferentially benefit TF_score-high patients. Restoring FOXO3, or IRF1 activity, through PI3K/AKT inhibition, epigenetic modulators (histone-deacetylase or LSD1 inhibitors), or interferon agonists—may concurrently lower stemness and enhance antigen presentation, priming tumors for checkpoint blockade (23). Finally, synthetic-lethal screening could uncover druggable co-factors unique to the MYC/E2F/FOXO3 axis, offering precision-therapy combinations guided by TF_score.
Several limitations merit acknowledgement. First, retrospective microarray data may lack dynamic range compared with modern RNA-seq; however, our external and single-cell validations mitigate this concern and underscore platform-independent robustness. Second, the OmniPath regulon, while comprehensive, may include false-positive edges; sensitivity analyses with DoRothEA and TRRUST regulons yielded concordant clusters, suggesting minimal bias. Third, immune-response predictions rely on in-silico surrogates; prospective correlative studies in patients receiving dinutuximab are needed to confirm TF_score as a predictive biomarker. Fourth, functional studies (CRISPR or degron-based knockdown of MYC/E2F, rescue of FOXO3/IRF1) were beyond the scope but will be essential to prove causality. Last, because both discovery (North American) and validation (Central European) cohorts were largely of European ancestry, the TF_score requires evaluation in ethnically diverse populations, including Asian and African cohorts, to confirm transferability of cut-off values and prognostic effect sizes.
Future work will focus on the following directions. (I) We will implement TF_score analysis prospectively in our institutional NB registry to confirm its prognostic value in a real-world setting. (II) Patient-derived xenograft models stratified by TF_score will be used to evaluate the efficacy of MYC/E2F-directed therapies and immune-checkpoint inhibitors, generating pre-clinical data for biomarker-guided trials. (III) A lightweight R/Shiny application will be released to allow diagnostic laboratories to compute TF_score directly from RNA-seq files, facilitating multi-center validation and routine adoption.
In summary, TF-activity profiling provides a practical TF_score that enhances current risk classification and offers a translational framework for personalized treatment in NB.
Conclusions
In summary, TF activity profiling allowed us to define biologically distinct NB subtypes and to derive a TF_score that robustly stratifies prognosis and can be implemented as a simple nomogram. Because TF_score can be calculated from routine transcriptomic data, it has the potential to complement existing risk stratification and to inform future translational and therapeutic work in NB.
Acknowledgments
None.
Footnote
Reporting Checklist: The authors have completed the TRIPOD reporting checklist. Available at https://tp.amegroups.com/article/view/10.21037/tp-2025-aw-739/rc
Peer Review File: Available at https://tp.amegroups.com/article/view/10.21037/tp-2025-aw-739/prf
Funding: This work was supported by
Conflicts of Interest: All authors have completed the ICMJE uniform disclosure form (available at https://tp.amegroups.com/article/view/10.21037/tp-2025-aw-739/coif). The authors have no conflicts of interest to declare.
Ethical Statement: The authors are accountable for all aspects of the work in ensuring that questions related to the accuracy or integrity of any part of the work are appropriately investigated and resolved. The study was conducted in accordance with the Declaration of Helsinki and its subsequent amendments. All datasets used in this study were obtained from publicly available repositories and contained only de-identified data; therefore, separate institutional review board approval and informed consent were not required.
Open Access Statement: This is an Open Access article distributed in accordance with the Creative Commons Attribution-NonCommercial-NoDerivs 4.0 International License (CC BY-NC-ND 4.0), which permits the non-commercial replication and distribution of the article with the strict proviso that no changes or edits are made and the original work is properly cited (including links to both the formal publication through the relevant DOI and the license). See: https://creativecommons.org/licenses/by-nc-nd/4.0/.
References
- Nong J, Su C, Li C, et al. Global, regional, and national epidemiology of childhood neuroblastoma (1990-2021): a statistical analysis of incidence, mortality, and DALYs. EClinicalMedicine 2025;79:102964. [Crossref] [PubMed]
- Monclair T, Brodeur GM, Ambros PF, et al. The International Neuroblastoma Risk Group (INRG) staging system: an INRG Task Force report. J Clin Oncol 2009;27:298-303. [Crossref] [PubMed]
- Bilke S, Chen QR, Wei JS, et al. Whole chromosome alterations predict survival in high-risk neuroblastoma without MYCN amplification. Clin Cancer Res 2008;14:5540-7. [Crossref] [PubMed]
- Chang J, Lin L, Zhang W, et al. Genetic variants of m(1)A modification genes and the risk of neuroblastoma: novel insights from a Chinese case-control study. Hum Genomics 2025;19:50. [Crossref] [PubMed]
- Guan Q, Lin H, Hua W, et al. Variant rs8400 enhances ALKBH5 expression through disrupting miR-186 binding and promotes neuroblastoma progression. Chin J Cancer Res 2023;35:140-62. [Crossref] [PubMed]
- Otte J, Dyberg C, Pepich A, et al. MYCN Function in Neuroblastoma Development. Front Oncol 2020;10:624079. [Crossref] [PubMed]
- van Limpt V, Schramm A, van Lakeman A, et al. The Phox2B homeobox gene is mutated in sporadic neuroblastomas. Oncogene 2004;23:9280-8. [Crossref] [PubMed]
- Shim J, Goldsmith KC. A New Player in Neuroblastoma: YAP and Its Role in the Neuroblastoma Microenvironment. Cancers (Basel) 2021;13:4650. [Crossref] [PubMed]
- Hoene V, Fischer M, Ivanova A, et al. GATA factors in human neuroblastoma: distinctive expression patterns in clinical subtypes. Br J Cancer 2009;101:1481-9. [Crossref] [PubMed]
- Badia-I-Mompel P. decoupleR: ensemble of computational methods to infer biological activities from omics data. Bioinform Adv 2022;2:vbac016. [Crossref] [PubMed]
- Wilkerson MD, Hayes DN. ConsensusClusterPlus: a class discovery tool with confidence assessments and item tracking. Bioinformatics 2010;26:1572-3. [Crossref] [PubMed]
- Zhang X, Shi M, Chen T, et al. Characterization of the Immune Cell Infiltration Landscape in Head and Neck Squamous Cell Carcinoma to Aid Immunotherapy. Mol Ther Nucleic Acids 2020;22:298-309. [Crossref] [PubMed]
- Yoshihara K, Shahmoradgoli M, Martínez E, et al. Inferring tumour purity and stromal and immune cell admixture from expression data. Nat Commun 2013;4:2612. [Crossref] [PubMed]
- Steen CB, Liu CL, Alizadeh AA, et al. Profiling Cell Type Abundance and Expression in Bulk Tissues with CIBERSORTx. Methods Mol Biol 2020;2117:135-57. [Crossref] [PubMed]
- Becht E, Giraldo NA, Lacroix L, et al. Estimating the population abundance of tissue-infiltrating immune and stromal cell populations using gene expression. Genome Biol 2016;17:218. [Crossref] [PubMed]
- Charoentong P, Finotello F, Angelova M, et al. Pan-cancer Immunogenomic Analyses Reveal Genotype-Immunophenotype Relationships and Predictors of Response to Checkpoint Blockade. Cell Rep 2017;18:248-62. [Crossref] [PubMed]
- Jiang P, Gu S, Pan D, et al. Signatures of T cell dysfunction and exclusion predict cancer immunotherapy response. Nat Med 2018;24:1550-8. [Crossref] [PubMed]
- Linstra R, Stappenbelt C, Bakker FJ, et al. MYC controls STING levels to downregulate inflammatory signaling in breast cancer cells upon DNA damage. J Biol Chem 2025;301:108560. [Crossref] [PubMed]
- Westwood JA, Matthews GM, Shortt J, et al. Combination anti-CD137 and anti-CD40 antibody therapy in murine myc-driven hematological cancers. Leuk Res 2014;38:948-54. [Crossref] [PubMed]
- Henssen A, Althoff K, Odersky A, et al. Targeting MYCN-Driven Transcription By BET-Bromodomain Inhibition. Clin Cancer Res 2016;22:2470-81. [Crossref] [PubMed]
- DuBois SG, Mosse YP, Fox E, et al. Phase II Trial of Alisertib in Combination with Irinotecan and Temozolomide for Patients with Relapsed or Refractory Neuroblastoma. Clin Cancer Res 2018;24:6142-9. [Crossref] [PubMed]
- Tian J, Wang J, Li S. Advances in the treatment of solid tumors in children and adolescents. Cancer Innov 2023;2:131-9. [Crossref] [PubMed]
- Mills CM, Turner J, Piña IC, et al. Synthesis and evaluation of small molecule inhibitors of LSD1 for use against MYCN-expressing neuroblastoma. Eur J Med Chem 2022;244:114818. [Crossref] [PubMed]

