henhaobio
Sign in

bioRxivCancer BiologyPreprint

Gene Regulatory Network Inference reveals tcf4 as a key a player in neuroblastoma gene expression circuitry

Catherine Koering, Elodie ValliniD, Franck PicardiD, Sandrine Gonin-GiraudiD, Olivier GandrilloniD

PDFDOICodeData
In-depth full-text digest available · Jump to digest ↓

Abstract

Neuroblastoma (NB), a pediatric cancer arising from disrupted sympathetic neuron differentiation, exhibits marked heterogeneity and limited therapeutic options. To better understand its molecular circuitry dynamics, we applied CardamomOT, a novel Gene Regulatory Network (GRN) inference framework, to single-cell RNA-seq data from patient-derived tumoroids. This approach models gene regulation via piecewise deterministic Markov processes, capturing transcriptional bursting and protein-mediated feedback, overcoming limitations of RNA velocity (e.g., gene independence and lack of biological time). We identified a continuous chromaffin-to-sympathoblast differentiation trajectory along which we selected 85 dynamically relevant genes enriched in cell cycle and DNA replication functions. Notably, 9 genes overlapped with those driving normal sympathoadrenal differentiation, underscoring tumor-normal tissue similarity. The inferred 85-genes network reproduced quite well experimental gene expression patterns in silico, and allowed to predict protein-level dynamics. Furthermore, it allowed to predict the effect of perturbations (both knock-out and overexpression) of hub genes (e.g., tcf4 and PLK1). We show that those perturbations significantly altered cell fate proportions in silico, with tcf4 KO increasing chromaffin-like cells and reducing proliferative late sympathoblasts. Predictions regarding tcf4 were tested using drug inhibition as a proxy for the gene KO. Using the BET inhibitor JQ1 indeed induced profound effect on the transcriptomic identity of our tumoroids. All of the 50 predicted tcf4 target genes were found to be significantly altered by JQ1 treatment. Finally cell fate proportions were also altered ex vivo closely resembling the predicted output. Our work therefore demonstrates that NB tumoroids retain a dynamic, differentiation-like architecture amenable to GRN modeling. Predicted druggable targets offer testable therapeutic avenues, including repurposing BET inhibitors or PLK1 inhibitors, potentially in combination.

In-depth summary

Claude Opus 5.5 · xhigh effort · 8 Nov 2026 · from the full text of version 2 (12 Oct 2026) · CC BY 4.0

Mechanistic stochastic gene regulatory network inference on neuroblastoma tumoroid scRNA-seq predicts TCF4 as a hub whose loss slows chromaffin-to-sympathoblast progression, supported by JQ1 treatment ex vivo.

Neuroblastoma, particularly high-risk disease, has few actionable targets, and earlier dynamic single-cell studies relied on RNA velocity, which fits genes independently and lacks an explicit biological time. Here an executable mechanistic gene regulatory network (GRN) fitted to patient-derived tumoroid data simulates differentiation and nominates perturbation targets such as TCF4 and PLK1.

Key results

  1. ShownSelecting genes by shifts in their expression distributions across pseudotime yielded 85 genes enriched for cell division and MCM genes, 9 shared with normal human sympathoadrenal differentiation; the inferred network shows hub-like outgoing connectivity, with TCF4 the most connected node.

    85 genes = intersection of 586 top Kantorovich-distance and 496 top delta-entropy genes; shared genes CHGA, DLK1, TH, HMGB2, IGFBP2, MDK, MEG3, RGS4, VIM; incoming degree near-Gaussian vs heavy-tailed outgoing degree; C4 feed-forward loops strongly enriched, C2 and I1 to I3 motifs underrepresented. · Results; Discussion; Figure legends

  2. SuggestedIn silico TCF4 knockout slows the chromaffin-to-sympathoblast progression, leaving many early sympathoblasts at the final time at the expense of late sympathoblasts; overexpression accelerates it, PLK1 and MKI67 act similarly but weaker, and TCF4 plus PLK1 double perturbation is additive.

    Perturbation of the 15 genes with most outgoing edges; WT vs perturbed composition compared by chi-square with Bonferroni correction (threshold 0.05); poorly connected CYGB and H1.0 had no effect; KO fixed protein at zero and basal parameter at -1000; simulated cells labelled by random forest (n_estimators=100). · Results; In silico perturbations; Figure legends

  3. ShownThe BET (bromodomain and extra-terminal) inhibitor JQ1, an indirect TCF4 proxy, reduced tumoroid cell numbers and remodeled transcriptomes beyond proliferation genes, shifting late sympathoblasts mainly toward early sympathoblasts, qualitatively matching predictions.

    10^-6 M JQ1: cell-number drop significant at 48 hours (n=5, Student's t test); ktest global p-value 4.10^-51 (2201 control vs 1532 treated cells); 120 induced and 484 repressed genes at adjusted p < 10^-3 and log fold change > 1, with induced-gene GO terms driven by histone genes; all 85 network genes affected; modest chromaffin gain. · Results; Tumoroids treatment; Figure legends

  4. SuggestedThe executable 85-gene model, simulated over 5 time points from 0 to 50 hours, approximates observed mRNA marginals, UMAP structure and cell-type proportions, and predicts protein trajectories splitting into rising and falling classes.

    Fit assessed by comparing simulated and experimental UMAPs and by computing distances between mRNA distributions; simulated repartition overestimated early sympathoblasts; networks from 3, 5 or 10 pseudotime slices all retained TCF4 and gave reasonable cell-type fits under Leiden or Louvain labels. · Results; Limitations of the study; Figure legends

  5. ShownTumoroid cells from one MYCN-amplified patient form a continuous trajectory from non-proliferating VEGFA+ chromaffin-like cells to MYCN+ proliferating sympathoblasts, with a side cluster mixing sympathoblast and Schwann cell precursor signatures.

    Phiclust clusterability index 0.9992 on 4345 cells; 8 Leiden clusters merged into chromaffin, early and late sympathoblast types; scVelo, CellRank and scFates agree on direction; Leiden vs Louvain adjusted Rand index 0.60 (clusters) and 0.66 (cell types); 90% of cells co-express PHOX2B and VIM. · Results; Discussion; Tumoroids generation and culture; Figure legends

Methods, data and reproducibility

Velocity pseudotime from scRNA-seq of patient-derived tumoroids was sliced into five pseudo-time points; 85 genes selected by Kantorovich distance and delta-entropy were fit with CardamomOT, a piecewise deterministic Markov model of bursty transcription with protein feedback, then simulated and perturbed in silico via HARISSA.

  • RevGel-seq scRNA-seq of neuroblastoma patient-derived tumoroids, two batches from one culture merged after QC · human (patient SJNBL013762; reads mapped to GRCh38) · n = 4345 cells after filtering (batches of 1744 and 2601 cells) · E-MTAB-16741
  • RevGel-seq scRNA-seq of tumoroids untreated or treated with 10^-6 M JQ1 for 24h · human · n = 2201 control cells; 1532 JQ1-treated cells · E-MTAB-16963
  • Orthotopic patient-derived xenograft cells from a 15 month-old boy (adrenal primary, bone marrow metastasis, MYCN amplification, ALK mutation) via St. Jude CSTN · human · n = one patient · SJNBL013762_X1
  • Subcutaneous PDX grown in female NSG-NOD mice to expand OPDX cells before tumoroid derivation · mouse (host) carrying human tumor cells
  • Tumoroid cell counts with or without JQ1 at 24 and 48 hours · human · n = n=5
  • CSTN Data Portal JQ1 dose-response data for neuroblastoma cells, including OPDX-derived primary cells and the SKNAS line · human · n = OPDX primary cells from 4 different patients
  • Previously published gene list from normal embryonic human adrenal gland sympathoadrenal differentiation, used for overlap · human
Methods & tools
CardamomOT (plos-cb-2026 version), HARISSA (stochastic GRN simulation), RevGel-seq (Scipio bioscience), STARsolo (Velocyto mode) with GRCh38, Seurat 5.4.0, scanpy, scVelo, CellRank and scFates, Phiclust, ktest (kernel-based differential testing), clusterProfiler, MultistageOT, scikit-learn RandomForestClassifier
Software & environment
CardamomOT plos-cb-2026 version; Seurat 5.4.0; STARsolo on GRCh38 (GRCh38.99.chr.gtf); inference settings SFT=10, stim=0.22, m=0.75, averaged over 10 parallel runs; simulation flags set in base.py with an inter_ref.csv reference GRN.
Assay conditions & compute
Inference run at the IN2P3 computing center as 10 parallel inferences; hardware, runtime and memory not reported.
Reproducibility gaps
Custom R scripts for motif enumeration and permutation tests have no stated repository; several checks are reported as not shown (batch UMAP, GeneCover/mixHVG gene selection, culture heterogeneity).
Limitations (6)
  • All results come from tumoroids of a single patient; the authors note that relevance to other patients is open given neuroblastoma heterogeneity, compounded by variability from PDX growth in different mice and tumoroid establishment. · stated by the authors
  • Gene selection relies on an arbitrary gene-number cutoff and ignores gene function, so it may omit relevant genes, misses fast events such as phosphorylation, and yields edges that may be indirect; lncRNA 'proteins' are abstract metavariables rather than real products. · stated by the authors
  • Time is imposed ad hoc by slicing velocity pseudotime at 12.5 hours per slice over a 50-hour process; the model overestimated early sympathoblasts, and whether motif enrichment reflects biology or inference bias remains unresolved. · stated by the authors
  • TCF4 knockout was not performed; JQ1 served as a proxy, and the authors note partial drug inhibition differs from full gene loss. The JQ1-experiment control contained very few chromaffin cells, reflecting difficulty standardizing ex vivo samples. · stated by the authors
  • JQ1 acts on BET bromodomains rather than on TCF4 directly, and it altered all 85 network genes, so the ex vivo shifts cannot be ascribed specifically to TCF4 without a genetic knockout or TCF4-target-restricted readouts. · noted in this summary
  • Internal inconsistencies: the abstract says TCF4 KO raises chromaffin-like cells whereas Results report early sympathoblast gains; motif tests cite 10^5 permutations in a legend but 1,000 Monte Carlo iterations in Methods; no test statistic is given for the JQ1 cell-type shift. · noted in this summary

Read this if you model gene regulatory dynamics from single-cell snapshots or study neuroblastoma differentiation states and want a mechanistic GRN workflow linking in silico perturbations to drug testing in tumoroids.

Citations
0

Reproducibility ledger

Nobody has reported trying to reproduce this paper yet. Did you replicate an experiment or run the code? Your report saves the next lab a week.

Sign in to file a reproducibility report.

Related papers

Similarity Map
SelectedSame fieldCross-field
THISGandrillon '26Szemes '18Guyer '22Banerjee '20Kildisiute '20Gómez-Muñoz '26
72% semantic matchSystems Biology
Olivier Gandrillon·F1000Research·Jul 23, 2026· 1 citations