Human tissue collection and multiome profiling
Patients with suspected GB provided preoperative consent at Cambridge University Hospitals NHS Foundation Trust, with written informed consent obtained under the 2000 Declaration of Helsinki and ethical approval from the local Cambridge Research Ethics Committee (REC 18/EE/0172). Tumors were resected under 5-aminolevulinic acid (5-ALA) guidance, sampled at multiple sites per tumor (targeting superior, anterior, posterior, and medial regions when possible), embedded in OCT, and snap frozen. Cryosections from each block were assessed for RNA integrity after extraction, and only those with an RNA integrity number greater than 7 (Agilent Tapestation) were serially sectioned for parallel mononuclear and spatial profiling. Nuclei were isolated from 50 µm cryosections by Dounce homogenization and Percoll gradient purification, then processed with the 10x Genomics Chromium Next GEM Single Cell Multiome ATAC + Gene Expression assay at 5,000 to 10,000 nuclei per reaction, two to three reactions per site. In total, 57 sites were treated for 12 tumors. Libraries were sequenced in pairs on a NovaSeq 6000, targeting a minimum of 100,000 read pairs per core per modality. Paired-end Visium libraries were prepared from adjacent sections, prestained with hematoxylin and eosin (H&E) as a reference, and sequenced on a NovaSeq 6000, targeting a minimum of 50,000 read pairs per point. Complete tissue processing and library preparation are described in our companion manuscript.14.
Multiome quality control, integration and annotation of cell state
Reads from the snRNA-seq and snATAC-seq libraries were aligned with Cell Ranger ARC to a custom reference genome, combining the 10x Genomics pre-mRNA reference GRCh38 and the Cell Ranger ARC ATAC genome (v.1.0.1). Nuclei were filtered based on the number of genes and unique molecular identifiers (UMIs) and mitochondrial fraction, and doublets were removed with Scrublet.69. ATAC data was processed with ArchR70 using thresholds on fragment counts and transcription start site (TSS) enrichment, then limited to barcodes passing RNA quality control, yielding 1,025,329 nuclei. Transcriptome counts were integrated into scVI30using patient, tumor site, reaction date, and cell cycle phase as batch covariates. Malignant populations were called from copy number profiles inferred with inferCNV71and clusters were annotated based on marker expression. We assigned malignant subgroups to states according to the genetic module notation according to Neftel et al.10projection on an atlas of first trimester brain development72 and enrichment of gene sets. A mapping between the Neftel States and ours is shown in Figure 1a of the extended data. The Atlas annotation has been finalized in the companion manuscript14 after the analyzes reported here were performed. Therefore, the cellular state assignments used here correspond to an earlier version, but the two annotations agree for 99.85% of malignant nuclei and for all neuronal and vascular microenvironment populations. The thematic activities by cell as well as the two annotations are available via Zenodo at https://doi.org/10.5281/zenodo.21873381 (ref. 73).
Epigenetic plasticity score
Chromatin accessible at a gene locus that is not yet transcribed indicates a primed configuration ready for potential activation. We used the degree of decoupling between chromatin accessibility and transcription at the marker gene level to quantify epigenetic plasticity, that is, the ability of a cell’s chromatin landscape to support transcriptional identities other than the one currently expressed. Following a previous study15we trained a classifier to predict cell states from gene expression profiles and applied it to gene accessibility profiles of the same cells. Cells whose accessibility is consistent with multiple transcriptional states are classified with low confidence and classified as more plastic.
Epigenetic plasticity scores were calculated on the metacell groupings of the patient atlas, constructed by Leiden clustering (Scanpy74tl.leiden(); resolution = 1) of the latent variable incorporating scVI RNA in each cell state, donor and batch stratum, resulting in a total of 16,423 metacells. Cell state-donor-batch groups containing fewer than 40 cells were directly considered a metacell. For each metacell, we calculated mean RNA counts and mean ArchR gene scores, the latter pooling ATAC signal across promoters, gene bodies, and distance-weighted distal elements. We normalized both matrices to library size and log transformed the RNA count matrix. A multinomial logistic regression model was trained on the expression of 3,160 highly variable genes to predict the 17 cellular states shown in Figure 1a (80% metacells as training fraction; 20% retained fraction; extended data, Figure 1b). The model was then applied to the scaled ATAC gene scores, and epigenetic plasticity scores were estimated as the Shannon entropy of the resulting class probabilities. Metacell scores were projected onto individual constituent cells for visualization (Fig. 1c). Variance components (Fig. 1f) were estimated by ordinary least squares, with cell state and donor as fixed effects, followed by ANOVA. Section 2 on Supplementary Methods provides additional information on factors contributing to variance in epigenetic plasticity and overall RNA-ATAC coupling.
scDORI
scDORI infers eGRNs from paired single-cell or single-nucleus RNA-seq and ATAC-seq data. Dimensionality reduction and eGRN network inference are typically considered as separate steps, and networks are inferred within discrete, predefined cell type annotations. scDORI unifies the two, such that the latent factors themselves have regulatory significance. An encoder maps each cell onto a low-dimensional representation, and mechanically constrained decoders reconstruct the observed accessibility and expression from this. These are the constraints that separate scDORI from conventional dimensionality reduction: a factor reconstructs data well only if it corresponds to a coherent program in which TFs bind to accessible elements to control nearby genes. We call these factors topics. Each subject is a signed enhancer-mediated network linking co-accessible peaks, co-expressed genes, and their predicted activating and repressive regulators, such that the targets and mode of action of a factor may differ between subjects depending on the accessible peaks. Each cell is modeled as a mixture of topics rather than being assigned to a cluster, so that networks vary continuously from cell to cell and a single topic can be active in multiple states. Training is done in mini-batches, which allows us to scale to millions of cells without aggregation. Detailed explanations of the scDORI components and their motivations, as well as the training schedule and empirical null links for filtering eGRN links, are provided in Supplementary Methods 1.1 to 1.9, with the rationale for the previous methods in Supplementary Methods 1.10 and 1.11. A schematic representation of the scDORI architecture is shown in Figure 2a of the extended data.
scDORI GB Patient Atlas Training
scDORI was trained on single-core ATAC-seq and RNA-seq profiles of 95% of patient atlas cells and evaluated on the retained 5%, using 182,677 highly variable peaks, 3,192 highly variable genes, and 195 TFs to learn 50 latent factors (subjects). During phase 2 training for GRN inference (module 4), we froze the encoder and modules 1–3 to retain the topic representations learned in phase 1 (see Extended Data, Fig. 2a and Supplementary Methods 1.6). Additionally, we used empirical subject-level TF expression for GRN pruning, calculated by averaging the TF expression across each subject’s cells as input to Module 4, rather than using the denoised TF estimates from Module 3. The full hyperparameters are provided in Supplementary Methods Table 3 and guidance on their selection can be found in Supplementary Methods 1.8. We verified that the subject programs inferred by scDORI were consistent across model restarts (three distinct random seeds) by calculating the Jaccard similarity between the top 500 genes of the most active subject in each GB state and that of baseline run 0 (extended data, Fig. 4d). Run 0 was used for all analyzes below. All subjects with peak activity greater than 0.1 in cells were considered for downstream analyzes (31 of 50). Subject activity across all states and cell types, as well as TF activity across all subjects, can be found in Supplementary Table 2. Dominant gene programs per subject were identified by gene set enrichment analysis (GSEA, uncoupler v1.6.0)75 by comparing the top 250 genes in the subject-gene matrix (Supplementary Methods 1.9.3) to all gene sets (size 15-500 genes) in subcategories H, C2, C5 and C8 of MSigDB (v.2023.1.Hs). The enriched pathways across subjects are summarized in Supplementary Table 2 and in the extended data, Figure 3a.
Comparative analysis of scDORI
We evaluated the performance of scDORI according to three complementary criteria: predictions of peak-gene links, the learned latent space and the ability of learned eGRNs to explain differential gene expression. Pic-gene links were compared to SCENIC+ on ENCODE76.77 cell lines IMR90, K562, and GM12878, using simulated single-cell multi-omics profiles derived from global RNA-seq and ATAC-seq data for these lines, as previously described18. Predictions were evaluated against the benchmark of candidate activator-gene interactions (BENGI)33a curated set of pic-gene interactions constructed from CRISPR disruptions, chromatin interaction assays, and expression quantitative trait locus (eQTL) data for these lines, with performance quantified by the area under the precision-recall curve (AUPRC). The learned latent space was compared to MOFA29 and SCENIC+18the latter represented by AUCell78 scores on the same set of TFs as scDORI, both run on metacell profiles due to memory demands. Subject scDORI activities were grouped at the metacell level for fair comparison. The biological relevance of each latent…
Gn Health
