Experimental method detailsHuman participants and ethical approval

Frozen GBM specimens that were diagnosed as ‘glioblastoma, IDH-wildtype’ according to the World Health Organization 2021 classification were collected at seven institutions, listed below. Collection and genomic profiling were approved by the institutional review board of each institution, and all patients provided informed consent accordingly. The tumors were collected from the MD Anderson Cancer Center, approved by the Institutional Review Board of MD Anderson Cancer Center under the protocol 2012-0441. Duke University cohort (n = 12)—the tumors were collected from the Preston Robert Tisch Brain Tumor Center Biorepository, Duke University Hospital. It was approved by the Institutional Review Board of Duke University under the protocol Pro0007434. Tokyo University cohort (n = 30)—the tumors were collected from the Department of Neurosurgery, Tokyo University Hospital. It was approved by the Institutional Review Board of Tokyo University Hospital under the protocol G10028. Pitié-Salpêtrière Hospital cohort (n = 28)—the tumors were collected from Pitié-Salpêtrière Hospital. It was approved by the Onconeurotek Brain Tumor Bank of the hospital Pitié-Salpêtrière certification 96-900. St. Michael’s Hospital cohort (n = 14)—the tumors were collected from the Division of Neurosurgery, St. Michael’s Hospital, Unity Health Toronto. It was approved by the Research Ethics Board of St. Michael’s Hospital, Unity Health Toronto, under the protocol REB 13-141, and all patients provided signed informed consent accordingly. Seoul National University (SNU) cohort (n = 20)—the tumors were collected from the Department of Neurosurgery, SNU Hospital. It was approved by the Institutional Review Board of SNU Hospital (approval H-2004-049-1116). NORLUX Neuro-oncology Laboratory (NORLUX) cohort (n = 13)—the tumors were collected from Centre Hospitalier de Luxembourg (CHL, Neurosurgical Department). It was approved by the National Committee for Ethics in Research, Luxembourg, under the protocol 201201/06. Cohorts were added to the Institutional Review Board protocol Dana–Farber/Harvard Cancer Center 10-417. Patients were male and female. Clinical information of the cohorts is summarized in Supplementary Table 1.

Statistics and reproducibility

No statistical method was used to determine the sample size. Data were excluded from select analyses when the number of malignant cells within a sample was low in at least one time point. Throughout the figures in this manuscript (Figs. 2–5 and Extended Data Figs. 1–3, 5 and 6), the boxplots reflect the following summary statistics: (1) the line splitting the box represents the median value; (2) the lower and upper edges correspond to the first and third quartiles (the 25th and 75th percentiles); (3) the upper whisker extends from the edge to the largest value no further than 1.5× interquartile range (IQR) from the edge; (4) the lower whisker extends from the edge to the smallest value at most 1.5× IQR of the edge; and (5) if the plots do not include jitter points representing the actual data, then data beyond the end of the whiskers are plotted individually as outlying points. The statistical analysis described in this work was done using R version 4.0.1 and above.

Nuclei isolation from frozen tissue

Protocol 1 (laboratory 1 and laboratory 2 for the cohorts from Duke University, MD Anderson Cancer Center, Tokyo University Hospital, St. Michael’s Hospital and Pitié-Salpêtrière Hospital)—nuclei from frozen tumor tissue were isolated as follows: tissue was thawed and mechanically dissociated in salt-Tris buffer (10 mM Tris–HCl (pH 7.5), 146 mM NaCl, 1 mM CaCl2 and 21 mM MgCl2) with 0.49% CHAPS (3-((3-cholamidopropyl) dimethylammonio)-1-propanesulfonate) (Millipore, 28300). Single-nuclei suspensions were filtered using a 40 μm strainer, centrifuged at 500g for 5 min and resuspended in salt-Tris buffer supplemented with 0.01% BSA (NEB, B9000S). Nuclear suspensions were inspected by microscope, counted using a hemocytometer and used for 10x Genomics workflow.

Protocol 2 (laboratory 3 for the cohorts from NORLUX and SNU cohorts)—tissue samples were thawed and mechanically dissociated in nuclei EZ lysis buffer (Millipore Sigma, NUC101) via Dounce homogenization. The solutions were incubated on ice for 5 min and mixed one to two times during incubation. Single-nuclei suspensions were filtered through a 70 μm strainer and centrifuged at 500g for 5 min at 4 °C, resuspended in nuclei EZ lysis buffer and incubated on ice for 5 min. The solutions were centrifuged at 500g for 5 min at 4 °C and resuspended in 1% BSA/0.2 U μl RNase inhibitor/PBS buffer (three times). For the final resuspension, the DAPI was added to the buffer, the solution was filtered through a 40 μm strainer, cells were counted on a Countess II automated cell counter (Thermo Fisher Scientific) and nuclei were taken into the 10x Genomics workflow.

10x Genomics for single-nucleus sequencing

The 10x Chromium Single Cell 3′ Reagent Kit v3 (10x Genomics, PN1000128) was used according to the manufacturer’s protocol. In brief, nuclei were loaded on the Chromium Chip (10x Genomics, PN12000120) with a target cell recovery of 6,000–8,000 nuclei and processed in the Chromium Controller. Single nuclei were partitioned into gel beads-in-emulsion (GEMs), followed by RNA reverse transcription with barcoding. Libraries were created by breaking GEMs and pooling barcoded fractions, cDNA amplification, fragmentation and attachment of sample index and adapter and sequenced on NextSeq500 or NovaSeq (Illumina). For NORLUX and SNU cohorts, nuclei were loaded on a Chromium chip with a target cell recovery of 6,000 nuclei for single-cell multiome ATAC and gene expression according to the manufacturer’s protocol. The gene expression chemistry for the 10x multiome possesses the same chemistry as the 10x v3 3′ snRNA-seq.

WES

WES for the samples from Duke University, MDACC, Tokyo University Hospital and St. Michael’s Hospital was performed as follows: DNA was extracted from each frozen tumor sample and blood sample corresponding to the patients using the DNeasy Blood & Tissue kit (Qiagen, 69504). The genomic DNA (100–250 ng) was acoustically sheared by an ultrasonicator (Covaris), targeting 150 bp fragments. Library preparation was performed by KAPA HyperPrep Kit (KAPA Biosystems, KK8504) followed by the enzymatic clean-up using AMPure XP beads (Beckman Coulter, A63882). Exome capture was performed using a custom exome bait (manufactured by Twist Biosciences to Broad Institute’s specification). Captured libraries were sequenced with 150-base-pair paired-end sequencing on a NovaSeq 6000 (Illumina). For the samples from Pitié-Salpêtrière Hospital, after the DNA was fragmented by a LE220 ultrasonicator (Covaris) and size selected, library preparation and capture were performed using the Twist Human Core Exome kit (Twist Bioscience) according to the manufacturer’s protocol. Sequencing was performed on a NovaSeq 6000 (Illumina) with 200-base-pair paired-end sequencing.

WGS

Newly generated whole-genome DNA-sequencing data were collected for a cohort of frozen samples from the NORLUX Neuro-oncology Laboratory. Briefly, DNA was extracted from each tumor sample using the AllPrep DNA/RNA Minikit (Qiagen, 80204) for samples with sufficient tumor tissue and matched normal blood when it was available. Note that the selected tissue for bulk DNA isolation was adjacent to the tissue used for single-nuclei data generation. Briefly, DNA was sheared to 400 bp using a LE220 ultrasonicator (Covaris) and size-selected using AMPure XP beads (Beckman Coulter, A63882). Whole-genome libraries were prepared and sequenced with 150-base-pair paired-end sequencing on a NovaSeq 6000 (Illumina). WGS data for the SNU cohort were prepared identically and were previously reported in ref. 9. Data for the SNU cohort are available on Synapse (https://www.synapse.org/glass).

Somatic variant detection and analysis

DNA-sequencing alignment, fingerprinting, somatic variant detection (Mutect2) and copy-number segmentation were performed in accordance with the Genome Analysis Toolkit (GATK) best practices using GATK 4.0.10.1, as previously described7,9. Briefly, whole-exome and whole-genome reads were aligned to the b37 genome (human_g1k_v37_decoy) using BWA MEM 0.7.17. DNA fingerprint analysis using ‘CrosscheckFingerprints’ (Picard) confirmed that all samples belonging to a patient came from the same individual, indicating that there were no sample mismatches in this study. Somatic mutations were detected using Mutect2 (v4.1.0.0) for tumor samples with a matched normal blood sample. A panel of normals was constructed for each sequencing batch to account for differences in DNA library prep (for example, whole genome versus whole exome) and used to filter out common artifactual and germline variants. Patient tumor samples without matched normal blood were analyzed for copy-number alterations and tumor-only SNV detection was performed with a panel of normal reference according to GATK best practices. Tumor-only somatic variants were used to identify hypermutant tumors. To identify samples with increased small deletion mutation burden, the following thresholds were applied: the recurrence-specific small deletion mutation burden needed to be greater than 0.2 Mut per Mb sequenced and greater than 0.1 Mut per Mb increase when comparing all small deletion variants in the recurrence versus the primary tumor. Mutational signature estimation was performed using COSMIC v3 signatures for the GLASS and CARE datasets according to previously published methods41. Briefly, all WGS samples were used to generate a mutational matrix and extract de novo signatures using SigProfilerExtractor (v.1.1.4)41 using maximum signatures set to 10 and NMF replicates set to 100. The de novo signatures were then deconvoluted to COSMIC v3.3 signatures. Finally, mutations for each sample were assigned to the most probable signature among the deconvoluted COSMIC signatures using the R package Palimpsest (v.2.0.0)42. Acquired SBS21 samples were denoted as those that were in the top quartile of SBS21 change at T2. For samples with matched blood, hypermutation was determined according to the cutoff of 10 Mut per Mb and further separated into de novo versus treatment-associated based on COSMIC mutational signatures. Treatment-associated hypermutation in samples with matched normal blood was classified by a longitudinal increase in mutation burden (that is, from 10 Mut per Mb) and SBS11 signal (alkylating agent-associated signature). For samples without matched blood, a single recurrent sample with a large increase in mutation burden at recurrence compared to the initial tumor (18.3-fold increase) and SBS11 signature was indicative of treatment-associated hypermutation. Genetic distance between time-separated samples from the same patient was defined as the bulk DNA proportion of private mutations (SNV) and differences in CNA segments between the two samples.

Single-cell/single-nucleus data processing of human glioma samples

For droplet-based snRNA-seq data, the Cellranger v3.1.0 pipeline provided by 10x Genomics was used for alignment (GRCh38 release 93) and to generate count matrices. We quantified the gene expression levels as \(\scriptstyle{E}_{i,j}=\log_2(\frac{{{\mathrm{CPM}}}_{i,j}}{10}+1)\), where CPMi,j refers to counts-per-million for gene i in sample j. We divided the TPM/CPM values by 10 as the size of single-cell libraries is estimated to be in the order of 100,000 transcripts, and we, therefore, would like to avoid inflating the expression levels by counting each transcript ~10 times. Following these normalization steps, initial filtering of low-quality cells was done based on low number of detected genes (less than 200 genes) or high expression of mitochondrially encoded genes (more than 3%). The majority of low-quality cells were removed downstream in case they could not be robustly classified as malignant or nonmalignant (see in later section on CNA analysis ‘Classification of cells into malignant or nonmalignant cell types’). Finally, we computed the average expression of each gene i as \(\scriptstyle\log_2((\frac{1}{n}{\sum }_{j=1}^{n}{{\mathrm{TPM|CPM}}}_{i,j}+1))\) and filtered out the lowly expressed genes by limiting the analyzed genes to the top 3,000 most highly expressed genes (using the Seurat v4.0.4 (ref. 43) function CreateSeuratObject). For the cells and genes that passed these quality control filters, we defined relative expression levels by centering the expression levels for each gene across all cells in the dataset as follows: \({{\mathrm{Er}}}_{i,j}={E}_{i,j}-\,\frac{1}{n}{\sum }_{k=1}^{n}{E}_{i,k}\), where n is the number of cells in the dataset.

Clustering and cell-type annotation

To facilitate clustering, cell-type annotation and scoring of droplet-based data, we used the R package Seurat (v4.0.4)43. Unique molecular identifier counts data were normalized and scaled (using the NormalizeData and ScaleData functions) and then clustered (by projecting the data to a lower dimensional space using principal component analysis and then running FindNeighbors and FindClusters functions). Initial classification of tumor-associated macrophages, T cells, oligodendrocytes and endothelial cells was based on the expression of known marker genes16,17,18,19,44, whereas the rest of the cells were annotated as presumed malignant. Following classification of cells as malignant or nonmalignant (see sections ‘Inferring CNAs from gene expression data’ and ‘Classification of cells into malignant or nonmalignant cell types’), the cells classified as nonmalignant were scored (using the AddModuleScore function) for known cell-type signatures and assigned with a cell type using the method explained in the section ‘Assignment of cells to states’. Finally, heterotypic doublets were filtered out using the R package scDblFinder45.

Inferring CNAs from gene expression data

CNAs were estimated as previously described46,47,48,49 using the function infercna from the R package infercna (an implementation of the original method presented in ref. 50, which is available at https://github.com/jlaffy/infercna) with default parameters. Briefly, the algorithm sorts the analyzed genes in each sample according to their chromosomal location and applies a moving average with a sliding window of 100 genes within each chromosome to the relative expression values. The scores computed for the cells classified as nonmalignant (oligodendrocytes, macrophages and endothelial cells) define the baseline of normal karyotype, and their average CNA values are used to center the values of all cells.

Classification of cells into malignant or nonmalignant cell types

In cells classified using single-nucleus droplet-based sequencing, we generally encountered lower CNA signal and continuous signal distributions related to (1) lower data quality relative to that observed in single-cell data and (2) inclusion of nonmalignant cell types such as ACs and neurons in the presumed malignant group, as these cell types could not be classified a priori as nonmalignant due to transcriptional similarity with the malignant cells.

We, therefore, used a multistep classification scheme for the droplet-based single-nucleus sequencing data.

Step 1: cell classification by detection of copy-number events

For each cell j, we computed the CNA signal on each chromosome separately as \({{\mathrm{CNA}}}_{j}^{C}=\frac{1}{n}{\sum }_{i=1}^{n}{{\mathrm{CNA}}}_{i,j}^{C}\), where n is the number of genes on chromosome C included in the analysis. We fitted a normal distribution (that is, a classifier) to the signal distributions of the reference cells for each chromosome (using the function fitdistr from the R package MASS v.7.3-56). This enabled assigning each cell with a P value for each chromosome using a two-sided z test (function pnorm from the R package stats v.4.1.1), reflecting whether the cell has a CNA signal that is substantially different than expected from nonmalignant cells. These P values were further corrected for multiple testing (for each chromosome) using the Benjamini–Hochberg (BH) method, and cells with BH-adjusted P value P value 

Step 2: exclude questionable cells based on CNA correlation

We define the CNA correlation as the Pearson correlation between each cell and the average CNA profile of the tumor sample the cell originated from. We limit the computation to the chromosomes on which ‘real’ events were detected and a control panel of ~1,000 genes from other chromosomes. Based on the distribution of CNA correlations across the cells classified as Malignant, we defined 0.5 as the lower threshold of malignancy. To determine the upper threshold for nonmalignancy, we performed SNV calling from the snRNA-seq data to estimate the algorithm misclassification rate. Briefly, for each sample with available WES data, we called SNVs for each cell using the Vartrix tool (https://github.com/10XGenomics/vartrix) from the FASTQ files. Misclassified cells were defined as cells harboring a malignant mutation that were classified as nonmalignant in step 1 of the algorithm. Following a cost-effectiveness analysis (misclassification rate versus percentage of excluded cells), we set 0.35 as the upper threshold for nonmalignancy. Finally, cells classified as ‘malignant’ with CNAcor  0.35 were reclassified as ‘Unresolved’ and excluded from downstream analysis.

Step 3: refine classification by setting confidence levels

To refine the classification, we defined confidence levels. Confidence level 1 includes cells classified as malignant that harbor a malignant SNV, cells in which at least 50% of the expected CNA events were detected, including chromosome 7 amplification or chromosome 10 deletion, and cells classified as nonmalignant with CNAcor 

In the analysis presented in this work, we used only cells in confidence levels 1 and 2.

Deriving metaprograms from gene expression data

To capture the heterogeneity between cells from the same cell type, we leveraged NMF51. NMF was performed on the relative expression values of each sample independently after setting negative values to zero. The NMF algorithm requires defining a priori the ‘k’ parameter, reflecting the expected number of latent features in the data. Because ‘k’ varies between samples and is largely unknown, we ran the NMF algorithm on each sample using different values (3–10), thereby generating 52 programs for each sample. Each of these NMF programs was summarized by the top 50 genes based on the NMF coefficients. Derivation of the metaprograms from the NMF programs is described thoroughly in ref. 24 and will be described here briefly. The method first filters out NMF programs that are not robust (do not recur within the sample it originated from or across several samples) or are redundant within the sample (that is, it highly overlaps with other NMF programs within that sample). Following this, the robust NMF programs are clustered according to the Jaccard similarity using a customized clustering algorithm that iteratively considers the degree of overlap between robust NMF programs and combines the highly overlapping ones to form a cluster. The top 50 most recurring genes in each cluster define a metaprogram.

Overall, the algorithm yielded 16 metaprograms that were annotated based on functional enrichment analysis (using Gene Ontology (GO) terms, mSigDB Hallmark gene sets and gene sets derived from normal brain development datasets). MP 16 included genes from various cell types, resulting in a metaprogram that fits almost all cells and thereby did not reflect any heterogeneity and was therefore excluded from analysis. This same procedure was independently repeated for each nonmalignant cell type.

Assignment of cells to states

Malignant cells were scored for the NMF metaprograms using the Seurat function AddModuleScore. To facilitate cell classification, we generated 20 shuffled expression matrices by sampling 5,000 cells each time and shuffling the expression values for each gene. We then scored each shuffled matrix for the NMF metaprograms, thereby yielding 100,000 normally distributed scores for each expression program. These served as null distributions for cell state classification. For each original cell, we computed a P value for each of the expression programs with a z test (R’s pnorm function) using the statistics of the null distributions that we previously generated. We adjusted all P values for multiple testing using the BH method. Each cell was classified into a specific state if the adjusted P value computed for that state was smaller than 0.05. Cells that achieved an adjusted P P 

Cells were assigned a ‘cycling’ state on top of their cellular state if they achieved an adjusted P 

Hybrids detection and classification

We define hybrids as cells scoring significantly (that is, with adjusted P 

To estimate the expected frequency of technical hybrids, we first computed for each pair of states, \(\left(A,{B}\right)\), the expected frequency of hybrids, defined as \({\mathrm{Exp}}\left(A,B\right)={\mathrm{Freq}}\left(A\right)\times {\mathrm{Freq}}\left(B\right)\), and the O/E ratio defined as \({\log }_{2}(\frac{{\mathrm{Obs}}(A,{B})}{{\mathrm{Exp}}(A,B)})\). As this is an overestimation of the expected frequency attributed to technical effects, we estimated this technical factor by averaging across O/E ratios of hybrid pairs detected less than expected by chance (that is, O/E O/tE ratio as \({\log }_{2}(\frac{{\mathrm{Obs}}(A,{B})}{{\mathrm{TF}}\times {\mathrm{Exp}}(A,B)})\). Finally, we filtered out insignificant hybrid pairs with large O/tE ratio (attributed to small numbers) using a one-sided Fisher test. The hybrid lineage model was generated using the R package igraph with cell states as nodes and edges connecting the states that form significant hybrid pairs. We did not include the ‘cycling’ state in this analysis because we don’t view ‘cycling’ as an independent state but rather as an additional feature that cells can have on top of their neurodevelopmental identity (that is, AC-like, OPC-like, etc.).

Baseline profile derivation

Each tumor sample was decomposed into seven pseudobulk profiles (one for each of the common states—AC-like, MES-like, hypoxia, GPC-like, OPC-like, NPC-like and NEU-like—given at least 25 cells classified to the particular state) by averaging across the normalized unique molecular identifier counts and log2-transforming. Genes were included if their mean log2-expression value across all pseudobulks was at least 1 and if their median variance (computed separately within state) was at least 2.5. Overall, 1,005 genes passed these filters. Following this, the pseudobulks were analyzed within the state to remove the state-specific signal using principal component analysis. We then derived six gene programs from each state (top and bottom 50 genes of the first three principal components), computed the Jaccard similarity index between each pair of programs and clustered the similarity matrix using hierarchical clustering. This revealed five distinct clusters from which we derived five consensus signatures by including genes that recurred in at least 25% of programs in each of the clusters (with a hard minimum of at least three times). By manual inspection and GO enrichment, we annotated the consensus signatures. Two signatures were excluded due to short length or high suspicion of reflecting a technical artifact. To exclude the possibility that ambient RNA has a role in driving the baseline profiles, we estimated per sample using the R package SoupX the contamination by ambient RNA. There was no difference in the estimated contamination level across the different baseline profiles as well as no significant overlap between the topmost contaminated genes per sample and the metaprogram or baseline profile signatures (data not shown due to lack of significance).

Assigning tumors to composition groups

To measure the differences in tumor composition, we generated for each tumor sample a composition vector reflecting the proportion of each cell type in the tumor sample. This enabled assigning tumors to the following three main composition groups: HMF (percentage of malignant cells >75%), IMF (percentage of malignant cells, 50–75%) and LMF (malignant cells 40%, or mixed in case of no dominant TME cell type). Similarly, we defined proportion vectors for the malignant cell states and assigned tumors to six groups according to the dominant cell state (at least 25% of cells assigned to the particular state)—AC, MES/hypoxia, GPC, OPC/NPC, neuron and mixed (in case of no dominant state). Tumors were classified according to the baseline profile with a maximal score and to the mixed category if the maximal score was less than 0.25.

Multilayer group definition

The association between each pair of features (A, B) was computed using a binomial test with the number of times the pair was observed as the number of successes, the number of tumor samples as the number of experiments and the frequency of feature A times the frequency of feature B as the expected probability of success (under the assumption that the features are independent). We then generated an undirected graph by defining each feature as a node and connecting the nodes with edges in case a statistically significant association (unadjusted for multiple testing) between the two features was observed.

Measuring transcriptional distance between matched pairs

To measure the transcriptional distance between matched pairs, we generated state-controlled pseudobulk profiles from each pair. State control was achieved by downsampling both samples to 25 cells from each state, given that both samples contain at least 25 cells from the particular state, resulting in samples that contain the same number of cells and are balanced across the different gene expression profiles of the cellular states. Following this, the pseudobulk profiles were computed by averaging for each gene across the selected cells, resulting in two pseudobulk profiles. Finally, the Euclidean distance between the two vectors was computed and served as the metric for transcriptional distance.

Survival analysis

Survival analysis was performed using the R packages survival v.3.3-1 and survminer v.0.4.9. Kaplan–Meier estimates were computed using the survfit function and plotted using the ggsurvplot function (survminer package). Statistical significance of the difference between the survival curves was computed using the log-rank test (implemented in the function survdiff). Cox regression analysis and forest plots were performed and plotted using the same packages (coxph and ggforest functions).

Differentially expressed genes between MGMT-low and MGMT-high pairs

To compute the differentially expressed genes between MGMT-low and MGMT-high pairs, we first computed the relative gene expression across time points in each group (that is, longitudinal expression difference), which is defined for each gene i as the average expression of the gene across all T2 samples in the group minus the average gene expression across all T1 samples in the group. Following, we computed the Euclidean distance between these two vectors. MGMT-low genes were defined as the top 100 genes with the largest distance between the groups and positive longitudinal expression difference in the MGMT-low group and negative longitudinal expression difference in the MGMT-high group. MGMT-high genes were defined as the top 100 genes with the largest distance between the groups and positive longitudinal expression difference in the MGMT-high group and negative longitudinal expression difference in the MGMT-low group.

Reporting summary

Further information on research design is available in the Nature Portfolio Reporting Summary linked to this article.