A single-cell based precision medicine approach using glioblastoma patient-specific models

Glioblastoma (GBM) is a highly lethal brain malignancy that remains refractory to standard-of-care (SOC) therapy, consisting of surgery, radiation (XRT), and temozolomide (TMZ), with median survival of 14–17 months. GBM tumors are complex ecosystems comprising normal cell types and diverse malignant tumor-cell states, and this intratumoral heterogeneity complicates treatment based on bulk genomic, epigenomic, and transcriptomic biomarkers that may not be uniformly present. Heterogeneity leads to tumor-cell populations with distinct phenotypes that respond differently to selective pressures; SOC treatment may select or induce drug-resistant states along poorly understood evolutionary trajectories. Because tumors are not routinely sampled immediately after SOC in clinical care, residual disease and the cell states driving recurrence are poorly characterized, necessitating tractable model systems that recapitulate baseline heterogeneity and treatment-emergent states. Precision medicine has expanded from mutation-based profiling to include bulk transcriptional profiling, but bulk analyses obscure intratumoral heterogeneity and often fail to identify master regulators driving cell-state biology that could be therapeutically exploited. The authors developed a framework to model and characterize non-genetic tumor-cell states across GBM progression, including disease onset, immediately post-SOC, and recurrence. Applying single-nucleus multi-omics (snRNA-seq, snATAC-seq) to an individual patient’s initial biopsy, untreated and SOC-treated patient-derived xenografts (PDX), and matched XRT-treated recurrent tumor at autopsy, they show that untreated and treated PDXs recapitulate phenotypic states observed in primary and recurrent tumors. Using systems biology and network inference (scSYGNAL/MINER), they identify mechanistic drivers of treatment-induced epigenetic state transitions and match potential therapeutics to active regulatory networks across disease stages. This provides a proof-of-concept for designing precision therapeutic regimens that account for intratumoral heterogeneity within a given patient.

Study design and samples: A single patient with IDH–wild-type WHO grade IV GBM underwent surgical resection with snap-frozen tissue processed for single-nucleus analysis. Fresh biopsy was orthotopically xenotransplanted into four NSG (NOD-SCID Il2rg−/−) mice to generate PDX ‘patient avatars’ split into control (n=2) and treatment (n=2) groups. Treated mice received SOC chemoradiation (TMZ 50 mg/kg IP daily ×5 days; X-ray irradiation at 115 cGy/min on same days). PDX tumors were harvested at 24 h (immediate post-treatment, IPT) and 72 h (late post-treatment, LPT) after treatment completion. The patient enrolled in a rapid autopsy program; recurrent tumor was collected (post-mortem interval ~8.75 h), snap frozen. This produced a longitudinal ‘time-lapse’ dataset: primary tumor, untreated PDX, SOC-treated PDX (24 h, 72 h), and matched recurrent tumor. Single-nucleus multi-omics: snRNA-seq and snATAC-seq were performed using 10x Genomics platforms. snRNA-seq used Chromium 3' v3 chemistry with targeted recovery of 12,000 nuclei and sequenced to ~50,000 reads/cell (28:8:0:91 cycles). snATAC-seq used SingleCell ATAC v1.1, targeting ~5,000 nuclei, sequenced to ~50,000 reads/cell (50:8:16:50). PDX conditions were hashed with TotalSeq-A oligo-tagged antibodies; demultiplexing used Seurat HTODemux. Doublets were identified with DoubletDecon (snRNA-seq; 50 runs consensus) for non-hashed datasets (primary, recurrent) and with ArchR’s synthetic doublet approach for snATAC-seq. Quality control and preprocessing: snRNA-seq raw data were processed with Cell Ranger 3.1 (GRCh38). QC filters per dataset included mitochondrial gene fraction ≤20%; counts/cell and genes/cell thresholds tailored per dataset (primary: counts 500–50,000; genes 500–10,000; PDX: counts 500–24,000; genes 500–7000; recurrent: counts 500–4000; genes 500–30,000). After QC: primary 4456 cells, PDX 4388, recurrent 350. Normalization used Seurat SCTransform with regression of mitochondrial effects. Dimensionality reduction used PCA (top 30 PCs), shared nearest neighbor (SNN) graphs, and Leiden clustering (resolution 0.8); visualization with UMAP. Cell type and tumor-cell annotation: CNS cell types were annotated using module scores of established markers (Seurat AddModuleScore). Tumor cells were identified by inferCNV-based copy-number signatures (Chr7 gain, Chr10 loss) and markers (SOX2, PTPRZ1, EGFR). Meta-module states (MES, AC, OPC, NPC) and MES–PN axis assignments followed Neftel et al. and Wang et al. frameworks. snATAC-seq processing: ArchR (v0.9.5) was used for preprocessing, QC, and analyses. QC thresholds: unique fragments >1000, TSS enrichment >4; fragment size periodicity and doublets assessed. Primary snATAC: 3407 cells post-QC (median 29,268 fragments/cell); recurrent: 1425 cells (median 8033). Dimensionality reduction used iterative latent semantic indexing (LSI) and SVD. GeneScore matrices were computed; label transfer to snATAC-seq used Seurat FindTransferAnchors to align to snRNA-seq references. Regulatory network inference: scSYGNAL (adapted SYGNAL) and MINER were applied to snRNA-seq tumor cells to infer regulons (co-regulated gene sets with TF/miRNA drivers) via biclustering and motif-enrichment, evaluate regulon significance by permutation (eigengene >95th percentile of null), compute per-cell regulon activity (eigengene values), and cluster regulons into transcriptional programs based on co-activity. Tumor cells were grouped into transcriptional network states by regulon/program activity profiles. Chromatin and TF activity: ArchR/chromVAR computed per-cell TF motif deviation scores and correlations to inferred gene expression (GeneScores) to identify TFs with differential motif activity. Consensus TFs were defined by overlap between scSYGNAL/MINER TF regulators and snATAC-seq motif activity evidence. Batch integration and projections: snRNA-seq tumor cells from primary and PDX were batch-integrated using Seurat (FindIntegrationAnchors, IntegrateData; CCA-based anchors). The integrated dataset included 7277 tumor cells with 6807 genes. Recurrent snRNA-seq profiles (6541 shared genes) were mean-centered/variance-normalized based on the integrated dataset and projected onto the integrated PCA/UMAP space; nearest neighbors to recurrent cells in UMAP were identified by Euclidean distance. Differential expression and enrichment: DEGs within SNN clusters were identified by Wilcoxon rank sum (|log2FC| ≥ log2(1.25), min.pct = 0.1, FDR ≤ 0.1). Functional enrichment used MSigDB Hallmark/C2/C4/C5/C6/C7 gene sets with clusterProfiler (FDR ≤ 0.1). Semantic similarity of regulons/programs used GO-based graph metrics (GOSemSim, Wang method; BMA aggregation). GSVA with permutation-based nulls assessed regulon enrichment in cells. Drug-target matching: Active regulons/programs across longitudinal stages (PT, IPT, LPT-A/B, REC) were matched to drugs targeting TFs or downstream genes using the Open Targets platform, focusing on agents in cancer clinical trials (Phases I–IV).

- Single-patient, longitudinal multi-omic model: Integrated single-nucleus RNA-seq (snRNA-seq) and ATAC-seq (snATAC-seq) from primary tumor, untreated and SOC-treated PDX avatars (24 h, 72 h post-therapy), and matched recurrent tumor at autopsy, enabling a ‘time-lapse’ view of GBM evolution. - Tumor-cell identification and heterogeneity: 4924 tumor cells (63.7% of total) identified in the primary tumor by Chr7 gain/Chr10 loss and GBM markers. Remaining cells represented expected microenvironmental types. - scSYGNAL/MINER regulatory landscape: Inferred 160 regulons (809 genes) co-regulated by ≥65 TFs and 141 miRNAs; 93 significant regulons comprising 52 non-redundant TFs regulating 454 targets. Of these TFs, 33 (63.4%, p=0.009) have prior GBM relevance, supporting biological validity. Regulons clustered into five transcriptional programs; tumor cells grouped into five transcriptional network states (SG-1..SG-5) with distinct biology (e.g., SG-4 hypoxia Pr-1 up, cell-cycle Pr-4 down; SG-3 proliferation up; SG-2 immune/cytokine; SG-1 angiogenesis; SG-5 mixed/transition with proliferation enrichment). - Multi-omic consensus TFs: snATAC-seq identified 141 TFs with differential motif activity; 103 showed motif deviation correlated with inferred expression. Overlap between scSYGNAL and snATAC defined seven consensus TFs (AR, TEAD1, RUNX1, RORA, EBF1, ZEB1, TCF4) implicated in YAP/Hippo, Wnt/β-catenin, EMT, and radiation resistance pathways. - PDX recapitulates primary phenotypes: Primary tumor formed 11 SNN clusters (3531 DEGs); PDX formed seven clusters (1731 DEGs) with enrichment of similar GBM-relevant processes. Nearly all primary clusters (except cl-6) were enriched for PDX cluster-specific DEGs. Cell-to-cell cosine similarity across transcriptional programs formed five mixed primary/PDX clusters with significantly higher PDX–primary similarities than primary–primary (p<0.05). Untreated PDX cells exhibited network states similar to 44.8% of primary tumor network states. - Longitudinal stages and recurrence mapping: Integrated UMAP revealed four stages: pre-treatment (PT), immediate post-treatment (IPT, 24 h), late post-treatment (LPT-A: 24/72 h; LPT-B: 72 h), and recurrent (REC). Recurrent tumor cells projected predominantly onto 72 h post-treatment PDX (PDX-tr2), with significant enrichment of PDX-tr2 regulons/programs in recurrent cells. - Dynamic regulon behavior patterns: Identified regulons with distinct longitudinal dynamics: • Selected-against (high early, reduced post-SOC): TEAD1-, TCF4-associated regulons. • Selected/induced (increased in REC): PRRX1 (TGF-β signaling) and IKZF2 regulons; IKZF2 potentially targetable via downstream CDK4 activity (CDK4/6 inhibitors). • Transient (peaking IPT/LPT): SP100 (regulon-7; enriched upregulated DEGs at 24 h only; includes PDGFC, MEOX2), and SOX5 (regulon-67) with druggable downstream nodes (TGF-β pathway, TTK). - Additional TFs from recurrent snATAC: RFX3 and RFX7 implicated in ciliogenesis and GBM resistance pathways; RFX7 not differentially expressed in snRNA-seq, highlighting multi-omic advantage. - Data integration scale: Batch-integrated dataset comprised 7277 tumor cells over 6807 genes; recurrent projection used 6541 shared genes. - Drug-matching exemplars: Potential therapeutics mapped to network vulnerabilities include galunisertib (TGF-β inhibitor) for PRRX1-related programs; CDK4/6 inhibitors (ribociclib, palbociclib) downstream of IKZF2; AR antagonists (enzalutamide, abiraterone) for AR-related programs; mineralocorticoid receptor antagonists (spironolactone, eplerenone) for NR3C2 (as depicted), and targets linked to SOX5/TTK.

This study presents a single-cell, multi-omic framework to reconstruct evolutionary trajectories of GBM within an individual patient across primary disease, immediate and late post-SOC treatment intervals, and recurrence. By integrating snRNA-seq and snATAC-seq with a systems biology network inference approach (scSYGNAL/MINER), the authors defined mechanistic transcriptional programs and network states that cut across conventional expression-based classifications, suggesting that targeting regulatory mechanisms could affect multiple phenotypic states simultaneously. The multi-omic integration provided consensus TF drivers not apparent from transcriptomics alone, underscoring the complementary nature of chromatin accessibility and gene expression in identifying regulatory determinants of cell-state transitions. PDX avatars recapitulated salient primary tumor phenotypes and, critically, captured transient and treatment-selected states immediately after SOC, a clinically unsampled window that likely seeds recurrence. Longitudinal categorization revealed dynamic regulon behavior—selected-against, induced/selected, and transient patterns—enabling rational mapping of candidate therapies to stage-specific vulnerabilities (e.g., TGF-β inhibition for PRRX1-associated induced states; CDK4/6 inhibition downstream of IKZF2; interventions against transient SP100/SOX5-governed programs during or immediately after SOC). The projection of recurrent cells onto 72 h post-treatment PDX further supports the model’s ability to anticipate recurrence-associated states. Overall, the findings address the research question by demonstrating that patient-derived, time-resolved, single-cell multi-omics combined with network inference can reveal actionable regulatory programs driving treatment-induced state transitions and recurrence, offering a precision medicine strategy tailored to intratumoral heterogeneity and temporal dynamics.

The study delivers a proof-of-concept precision medicine framework that: (1) reconstructs GBM tumor-cell evolutionary trajectories at single-cell resolution using patient-specific primary, PDX (untreated and SOC-treated), and recurrent samples; (2) infers mechanistic regulatory programs and network states via scSYGNAL/MINER and validates consensus TF drivers through snATAC-seq; (3) demonstrates that PDX avatars can model primary phenotypic states and SOC-induced transitions, including transient residual states; and (4) matches stage-specific network vulnerabilities to candidate therapeutics. Future work should functionally validate drug–target predictions, maintain and benchmark serial PDXs to ensure faithfulness to patient tumors, scale to larger cohorts to assess generalizability, and explore approaches to modulate non-druggable regulators (e.g., CRISPR/antisense) to translate regulatory network insights into clinical interventions.

- Single-patient case study limits generalizability; larger cohorts are needed to assess robustness across diverse GBMs. - Predicted drug–target pairings were not experimentally validated in this study due to limited biopsy material; functional testing in patient-specific models is required. - Lack of continuous, serially transplanted PDX maintenance precluded longitudinal in vivo validation; faithfulness of serial PDXs to original tumors requires systematic evaluation. - Many inferred regulators may not be directly druggable; validation and targeting may require genetic tools (e.g., CRISPR-Cas9, antisense RNA). - Multi-omic inference integrates orthogonal signals but still infers rather than directly measures TF activity; discrepancies between RNA and ATAC modalities can complicate driver identification.