Genomic Sequencing and Computational Analysis of IDH-Mutant Gliomas
DNA Sequencing: WGS and WXS Methods
Whole-genome sequencing (WGS) and whole-exome sequencing (WXS) data were newly generated for samples from the MD Anderson Cancer Center, Duke University Hospital, Henry Ford Hospital, Case Western Reserve University, Chinese University of Hong Kong, Mayo Clinic, Northwestern Medicine, Olivia Newton-John Cancer Research Institute, University of California, San Francisco, Saint Joseph’s Hospital and Medical Center, Seoul National University and Luxembourg Institute of Health at the Jackson Laboratory for Genome Medicine sequencing core.
All samples underwent 150-base-pair (bp) paired-end sequencing on an Illumina NovaSeq platform according to the manufacturer’s protocol. Previously published longitudinal data from GLASS5,12 and TCGA7 were also used; sequencing information for those datasets is available in the corresponding publications.
In total, sequencing data were contributed by 21 institutions: Columbia University, University of Florida, German Cancer Research Center, Pitié-Salpêtrière, Samsung Medical Center, Biobank Japan, University of South Wales, University of Sao Paulo and the 13 institutions listed above, including the Jackson Laboratory.
For the Seoul National University cohort, DNA and RNA were extracted concurrently from each tumour sample at the Jackson Laboratory for Genomic Medicine using QIAGEN’s AllPrep DNA/RNA Mini Kit (80204). DNA from blood samples was extracted using QIAGEN’s QIAamp DNA Mini and Blood Mini Kit (51104).
A total of 200 ng of DNA was fragmented to 400 bp using an LE220 focused-ultrasonicator (Covaris) and size-selected with Ampure XP beads (Beckman Coulter). DNA fragments underwent end repair, A-tailing and ligation of Illumina unique adapters using the KAPA Hyper Prep Kit from Roche (7962363001). The resulting whole-genome libraries were sequenced with 150-bp paired-end reads on the Illumina NovaSeq 6000.
Average sequencing coverage was 27× for WGS longitudinal samples and more than 100× for WXS longitudinal samples. Driver mutations were consistently supported by sufficient depth, ranging from 30×–80× for WGS and 80×–300× for WXS, in both primary and recurrent tumours. This minimized the effect of coverage variability on mutation detection.
Reference Cohorts: TCGA and AACR GENIE
GLASS data were compared with two reference cohorts: The Cancer Genome Atlas (TCGA)7 and AACR GENIE (version 18.5)25. Baseline clinical data for patients in the GLASS cohort were compared with data from TCGA patients with primary IDH-mutant gliomas, including age at onset and overall survival5,7,12.
Survival data for patients with oligodendrogliomas and astrocytomas were compared between the GLASS and GENIE cohorts. The GENIE cohort primarily contains panel-sequencing data and generally lacks copy-number data or coverage of TERT promoter hotspots.
Oligodendrogliomas in the GENIE cohort were therefore defined as IDH-mutant gliomas without alterations in TP53 or ATRX, based on high-coverage panel sequencing. IDH-mutant gliomas with alterations in these genes were classified as astrocytomas.
GLASS Sample and Patient Identifiers
A GLASS barcode system based on the TCGA barcode design was used to identify samples and their corresponding patients. The 24-character barcode contains information about the project, collection centre, patient identification, sample type, analyte type and sequencing platform. Detailed information about the barcode structure is available in a previous GLASS publication12.
Computational Pipelines for Genomic Analysis
GLASS pipelines for variant and copy-number calling were developed using Snakemake (version 5.2.2)46. The Genome Analysis Toolkit 4 (GATK4, version 4.1.0.0)47 suite was used unless otherwise specified. All data were analysed using homogeneous pipelines capable of processing raw FASTQ files.
Sequence Alignment, Quality Control and Variant Calling
Sequence alignment and preprocessing followed GATK Best Practices. Various tools were used to manage and sanitize binary alignment map files, assign and manage read groups, perform quality control, mark sequencing adapters, align sequences to the b37/hg19 genome (human_g1k_v37_decoy), merge and mark duplicates and perform base recalibration.
For variant detection, GATK4 was used to call germline variants from control samples and somatic variants from tumour samples and patient cohorts. Quality-control measures included assessment of cross-sample contamination and read-orientation artefacts.
During variant post-processing, BCFTools (version 1.9)48 normalized, sorted and indexed variants to generate a consensus variant call format (VCF) file. This file was annotated and used to compare allele frequencies between related samples. Single-sample variants were overlaid on multi-sample calls to infer individual sample mutations. Further details are available in the previously published GLASS paper5.
Copy-Number Alteration Calling
Copy-number identification followed GATK Best Practices, using different methods for WGS and WXS data. For WGS data, the genome was divided into 10-kb bins. For WXS data, overlapping regions from different capture kits were merged.
Both approaches included read counting within bins, subsetting germline variants and separating the cohort by sequencing centre to create a panel of normals. Read counts were denoised, intervals were annotated for GC content and results were inspected for quality.
Copy-number calling used a method based on GATK CallCopyRatioSegments and GISTIC (version 2.0.23)26. Neutral segments were identified, and amplifications or deletions were classified using weighted means and standard deviations. Chromosome arms were analysed to establish high-level thresholds.
Gene-level copy numbers were called by intersecting gene boundaries with segment intervals and comparing the weighted non-log2 copy ratio with these thresholds. These procedures have also been described in detail previously5.
Molecular Cohort Selection and Quality Control
Stringent quality-control measures were applied to retain time-separated surgeries for molecular analyses (Extended Data Fig. 1b). The criteria were:
QC1: A tumour sample was available and WGS and/or WXS data had been generated from the sample.
QC2: Fingerprinting confirmed sample identity using the Picard CrosscheckFingerprints tool49.
QC3: Sequencing coverage was at least 8× for WGS samples and at least 28× for WXS samples. Mutation frequency could not be zero or abnormally elevated, which would indicate possible FFPE artefacts.
QC4: At least 3 months separated the initial and recurrent surgeries.
QC5: Patients with only one surgery remaining were excluded. The analysis was restricted to two distinct surgeries per patient, using the first and last available surgery.
QC6: Copy-number profiles had to meet quality thresholds without excessive noise, unexpected genomic stability or insufficient signal. These profiles were independently assessed and jointly reviewed by three authors (T.C., K.C.J. and F.S.V.).
Single-Nucleotide Variant and dN/dS Analysis
dN/dS ratios were estimated using dNdScv (version 0.1.0)50, with the default and recommended parameters applied to all mutations in initial and recurrent tumour samples. Analyses were also performed for each mutational fraction: mutations unique to initial samples, mutations unique to recurrent samples and shared mutations.
These analyses were performed separately for the two types of IDH-mutant glioma. Adequate sequencing depth and reciprocal mutation calling between paired samples minimized coverage-related false negatives. Because the analysis compared shared mutations with sample-specific mutations, coverage variation was unlikely to affect the qualitative patterns observed.
Tumour Clonality Analysis
PyClone-VI51 was used to assess the clonal status of all single-nucleotide variants in the dataset. Variants with a minimum coverage of 10× in at least one sample from each paired set were included.
PyClone-VI was run in beta-binomial mode using default parameters. Mutations with a cancer-cell fraction greater than 0.9 were classified as clonal, whereas mutations below this threshold were considered subclonal.
GISTIC Analysis of Significant Copy-Number Changes
GISTIC26 was run on the public GenePattern server using segmentation data from primary and recurrent tumours. Each type of IDH-mutant glioma was analysed separately to identify significant broad and focal copy-number changes.
The analysis used the GeneGISTIC method, with broad-event detection and gene-collapse values set to extreme. Additional parameters included a focal-length cut-off of 0.5 or higher, a confidence level of 99% and a chromosome-arm adjustment capped at 1.0. Copy-number driver genes were identified from the GeneGISTIC results.
Aneuploidy and Fraction of Genome Altered
Aneuploidy was calculated at three levels for all samples: chromosome-arm level, segment level and genome-wide level.
For arm-level aneuploidy, a previously described method was used52. Adjacent segments with identical arm-level calls (−1, 0 or 1) were merged into a single segment. The proportion of each chromosome arm covered by the merged segment was calculated. Segments covering more than 80% of an arm were classified as −1 (loss), 0 (neutral), +1 (gain) or NA if no segment covered at least 80% of the arm.
For each sample, the number of chromosome arms with non-neutral events was counted. As a secondary method, average log2(copy-number ratios) were calculated to compare arm-level aneuploidies between groups. Thresholds of ±0.1 were used for significant arm-level gains or losses at any time point, and thresholds of ±0.05 were used for significant arm-level changes, including delta gains or losses.
For segment-level calculations, each segment was assigned a gain, loss or neutral status according to a log2 copy-ratio threshold of ±0.2. Segments with an absolute log2 ratio greater than 0.2 were considered altered. The total number of altered segments, including gains and losses, was counted for each sample.
Thresholds for high-level amplification and homozygous deletion were determined for each sample from segmented copy-number data. The weighted mean and standard deviation of non-log2 copy ratios among neutral segments (0.9–1.1) were calculated, with segment length used as the weighting factor. Iterative filtering was performed twice to remove outliers. Amplification and homozygous deletion thresholds were then defined as the twice-filtered weighted mean ± 2 standard deviations.
At the genome-wide level, the fraction of the genome altered (FGA) was calculated using the formula: FGA = 1 – (total base pairs in neutral segments/total callable base pairs).
For WGS samples, the callable genome size was assumed to be 3 × 109 bp. For WXS samples, callable size was set to 5 × 107 bp based on the exome-capture target region. FGA values were calculated for each sample using platform-specific copy-number profiles and normalized to a 0–1 scale. WGS and WXS represented mostly non-overlapping tumour sets.
Structural Variant Detection from WGS Data
Structural variants (SVs) were called only from WGS data using SvABA (version 1.2.0)53 against the hs37d5 reference genome. The analysis incorporated dbSNP build 138 (b37) and excluded ENCODE hg19 blacklist regions.
Somatic SVs were extracted from SvABA VCF files and parsed in R. Break-end records were resolved to obtain partner coordinates and strand orientation. Only autosomes and chromosome X were retained; events on decoy, random or alternative contigs were excluded.
Reciprocal break ends were paired to reconstruct complete events. Events were classified as deletions, duplications, head-to-head inversions, tail-to-tail inversions or translocations. Duplicate mates were collapsed into a single representative event per sample.
Chromothripsis Detection
Chromothripsis events were identified using ShatterSeek (version 1.1)54. Copy-number segments were derived from GATK segmentation output for WGS aliquots, and somatic SVs were called using SvABA. Only samples with both SV and copy-number data were retained for chromothripsis analyses.
Copy-number segment log2 ratios were converted to integer total copy numbers and capped at 20. Adjacent segments with identical copy numbers were merged.
Chromothripsis was assessed using break-point clustering, oscillating two-state total copy-number patterns and fragment-join randomness. Chromosome-level summaries were used to define confidence categories.
High-confidence events required more than five intrachromosomal SVs, more than six oscillating segments, significant break-point enrichment or clustering (P < 0.05) and non-random fragment joins (P > 0.05). Low-confidence events met similar criteria, with four to six oscillating segments. Alternatively, high-confidence events required more than two intrachromosomal SVs, more than three interchromosomal events and more than six oscillating segments.
Kataegis Analysis
Kataegis events were called using SigProfilerSimulator and SigProfilerClusters (version 1.2.2)55, as previously described.
For each sample, 100 simulated genomes were generated under a pentanucleotide context model while preserving chromosome-level mutation burdens and transcriptional strand bias. These simulations were used to derive empirical intramutational distance thresholds.
Cluster detection incorporated genome-wide correction for regional mutation density using 1-Mb windows. Variant allele frequencies were considered during subclassification.
Kataegis results were selected from the outputs, and mutations sharing a cluster-group identifier were aggregated and defined as a single kataegis event. The number of events per sample was counted to establish kataegis burden as log2(kataegis count + 1).
Kataegis analyses were restricted to non-hypermutator samples to avoid violating genome-wide background mutation-rate assumptions.
Whole-Genome Duplication Analysis
ASCAT was used in multi-sample mode to derive patient- and sample-level whole-genome duplication (WGD) status56. Denoised copy ratios and tumour and normal allelic count files from the GATK CNV caller were transformed into tumour and normal B-allele frequency and tumour and normal logR files for ASCAT input.
ASCAT was run in case-based multi-sample mode. WGD was called from the ASCAT solution as part of the reported summary metrics.
To avoid overestimating WGD, a second filter was applied. Samples predicted by ASCAT to have undergone WGD were retained only when the major copy-number fraction assessed by Sequenza was at least 0.5.
Copy-Number Signature Analysis
Allele-specific copy-number profiles from ASCAT were used as input for copy-number signature analysis with the SigProfilerToolkit55.
Sample-level 48-feature copy-number matrices were generated using SigProfilerMatrixGenerator. Copy-number signatures were assigned using SigProfilerAssignment55 and the COSMIC version 3.5 reference copy-number signature set38. Absolute and relative signature activities were exported for downstream analyses.
Extrachromosomal DNA and Focal Amplifications
Complex focal amplifications and extrachromosomal DNA (ecDNA) were detected using AmpliconSuite57, which integrates AmpliconArchitect and AmpliconClassifier.
Case-based multi-sample mode was used with default AmpliconSuite parameters. The output categorized amplicons into predefined classes, including ecDNA, breakage–fusion–bridge cycles and linear amplification. Samples were classified as either lacking or containing ecDNA.
Mutational Signature Analysis
Mutational signature analysis was performed in three steps using previously published methods28,58.
First, WGS data from all samples were used to create a mutational matrix and extract de novo signatures with SigProfilerExtractor (version 1.2.1)58. The analysis used the b37/hg19 reference genome, a maximum of 10 signatures and 100 non-negative matrix factorization replicates.
Seven SBS and seven ID de novo signatures were extracted from the oligodendroglioma cohort, while ten SBS and six ID de novo signatures were extracted from the astrocytoma cohort. All had average stability greater than 0.8.
Extracted signatures with a cosine similarity of at least 0.8 to a reference signature were deconvoluted. Signatures with a similarity below 0.8 were retained as novel signatures (Supplementary Table 12).
Second, de novo signatures were deconvoluted to a curated list of reference signatures for SBS, DBS and ID. The reference list combined two public mutational-signature databases: the COSMIC database (https://cancer.sanger.ac.uk/signatures)28 for SBS signatures 1–96 and all ID signatures, and the Signal database (https://signal.mutationalsignatures.com/explore/main/cancer/signatures)29 for SBS signatures 97 and higher.
To ensure confidence in signature-assignment probabilities, only WGS data from the cohort were used, consistent with the WGS-based curated reference.
Finally, the deconvoluted reference signatures identified in the second step were used to estimate mutational-signature exposure in each sample. The R package Palimpsest (version 2.0.0)59 was then used to assign the most likely signature to each mutation in each sample.
Four cases were excluded from downstream mutational-signature analyses because the initial sample, but not the recurrent sample, had a mutation burden exceeding 10 mutations per Mb.
Tumour Mutational Burden and Hypermutation
Mutational burden was defined as the number of mutations per Mb sequenced, requiring a minimum base coverage of 15×. Mutational burden was calculated for each sample and displayed on a common scale from 0–10 for non-hypermutated tumours.
WGS and WXS represented distinct tumour sets rather than matched samples sequenced on both platforms.
DNA hypermutation was defined in recurrent tumours by both more than 10 mutations per Mb and more than 25% of recurrence-specific mutations attributable to the COSMIC mutational signature SBS11.
Tumours with fewer than 10 mutations per Mb were classified as non-hypermutant. Hypermutation status was left unassigned for tumours with more than 10 mutations per Mb but fewer than 25% of recurrence-specific mutations attributable to SBS11.
Mutual Exclusivity of Genomic Alterations
Mutual exclusivity of genomic alterations was assessed using the DISCOVER algorithm (version 3.1.7)21. DISCOVER is a discrete independence statistic that controls for observations with varying event rates.
The algorithm calculated tumour-specific probabilities of gene alterations while preserving the observed marginal totals for each gene and tumour. Expected co-occurrence rates between gene pairs were derived using a Poisson binomial model that accounted for these probabilities.
Significantly fewer co-alterations than expected under the null model were interpreted as evidence of mutual exclusivity.
Assessment of Sequencing-Centre Batch Effects
Potential sequencing-centre-related batch effects were assessed by testing associations between treatment centres and recurrent single-nucleotide variants, copy-number alterations and mutational signatures.
Tumour-type- and hypermutation-status-stratified logistic regression models showed no significant centre-specific effects.
Overall Survival and Recurrence Analysis
Survival end points included overall survival from diagnosis and time to recurrence. Overall survival from diagnosis was defined as the interval between the date of primary surgical resection and the date of death, or the date of last follow-up if censored.
Time to recurrence was defined as the interval between the date of primary surgical resection and the date of first recurrence surgery. No censoring was required for time to recurrence because all patients had undergone a secondary surgery by study design.
Overall survival from birth was calculated by adding age at diagnosis to overall survival from diagnosis. Post-recurrence survival was defined as the interval from the date of the quality-control-passed last recurrence surgery to death or censoring at last follow-up.
Kaplan–Meier survival curves were generated and compared using the log-rank test with the survival and survminer R packages. Cox proportional hazards regression was used to estimate hazard ratios. Variables with a likelihood-ratio-test P ≤ 0.10 were considered for inclusion in multivariable models.
Schoenfeld residuals were examined to verify the proportional-hazards assumption. Variables violating this assumption were excluded. Prognostic models were refined through stepwise backward elimination, and model performance was quantified using Harrell’s concordance index with 1,000 bootstrap resamples.
Statistical Methods and Significance Testing
All statistical analyses were performed using R (version 4.3.2) and Python (version 3.10.13).
For cohort-level comparisons that did not account for patient-level pairing, unpaired categorical variables were assessed using Fisher’s exact test or the χ2 test. Unpaired continuous variables were assessed using the Wilcoxon rank-sum or Kruskal–Wallis test, as appropriate.
For cohort-level analyses evaluating patient-level discordance over time, paired categorical variables were compared using McNemar’s test and paired continuous variables using the Wilcoxon signed-rank test.
Statistical significance was defined as two-tailed P < 0.05 unless otherwise specified.
Ethics Statement and Informed Consent
The GLASS Consortium is an international observational study designed to investigate the molecular response to treatment in adults with glioma. Data were pooled from different sources.
The protocols used to collect and sequence specimens, and to collect and analyse patient data, were approved by the Yale University institutional review board (2000036542, 2000041140, 2000039350, 2000036441, 2000035189, 2000034831, 2000035190, 2000035426 and 2000035453) and as previously reported5,12.
Anonymized patient and sample identifiers were generated as previously reported5. All participants provided written informed consent before enrolment. The study was conducted in accordance with the Declaration of Helsinki and applicable institutional and national ethical guidelines.
Reporting Summary
Further information about the research design is available in the Nature Portfolio Reporting Summary linked to this article.
Source: www.nature.com