We investigated 88 retrospective leukemic bone marrow (BM) samples from patients diagnosed with pediatric B-cell ALL, treated according to NOPHO trials protocol at Uppsala University Hospital and Karolinska University Hospital in Sweden. The cases were selected to represent the mandatory genetic subgroups ( Supplementary Table 1 ) defined in the current protocol (n = 49), and the remaining samples consisted of B-cell ALL cases lacking recognized stratifying aberrations, i.e., the B-other group (n = 36), as well as three patients with Down syndrome ALL (DS-ALL). The blast count was above 50% in all but two samples: P111 with an ETV6::RUNX1 fusion present in 37% of the cells and P120 with a KMT2A rearrangement (KMT2A-r) present in 14% of the cells. The distribution of B-cell ALL samples across genetic subgroups is summarized in Table 1 . The samples were divided into an exploratory set (n = 58) and a validation set (n = 30). BM samples taken during follow-up, when the patients were in remission, were used as the source of germline DNA (n = 64) for the paired analysis.

At diagnosis, all samples were genetically characterized according to SoC as specified in the NOPHO protocols, which included chromosome banding analysis, fluorescence in situ hybridization (FISH) analysis to investigate the presence of t(12;21), t(9;22), KMT2A-r, intrachromosomal amplification of chromosome 21 (iAMP (21)) (Abbott, Vysis, Abbott Park, IL, USA), t(1;19) (MetaSystems Probes, Heidelberg, Germany), and dic(9;20) (Kreatech Diagnostics, Amsterdam, Netherlands) as well as array comparative genomic hybridization/single-nucleotide polymorphism (CGH/SNP) array analysis. For the purpose of the study, the genetic subgroup for the samples from the NOPHO 92 and 2000 trials (n = 22) was updated to SoC in the NOPHO 2008 trial, as the older trials investigated fewer aberrations (22, 23). Also, three of the B-other samples that harbored retrospectively detected ABL-class rearrangements were revised, as this subgroup was not investigated in the NOPHO 2008 trial.

BM samples were collected in 5-ml EDTA tubes, and genomic DNA (gDNA) was isolated using a Tissue kit an EZ1™ automated instrument (Qiagen, Hilden, Germany) or extracted manually from frozen cell pellets with the AllPrep DNA/RNA Mini Kit Qiagen (Qiagen, Hilden, Germany). DNA was stored at −20°C until use. The DNA from frozen BM samples taken at remission was extracted with a Gentra Puregene Blood core kit (Qiagen, Hilden, Germany).

For the samples in the exploratory set (n = 58), library preparation and sequencing were performed at the National Genomics Infrastructure, Science for Life Laboratory, Stockholm (n = 36) and Uppsala (n = 22), Sweden. For samples with sufficient amounts of DNA available (36/58), the libraries were prepared using 1 µg of DNA with the TruSeq DNA PCR-free (Illumina, San Diego, CA, USA) protocol, 350-bp insert size, and sequenced using the HiSeq X platform (Illumina), 2 × 150-bp paired-end to ~90× coverage for leukemia samples, and ~30× coverage for the paired germline sample. For samples where little DNA was available (n = 22), 100 ng of input DNA was used for library preparation with TruSeq DNA Nano protocol (Illumina).

The samples in the validation set (n = 30) were processed at Clinical Genomics, SciLifeLab, Stockholm. Libraries were prepared using the NxSeq^® AmpFREE Low DNA Library Kit (Lucigen, Biosearch Technologies, Petaluma, CA, USA) with 200 ng of DNA input and thereafter sequenced on the NovaSeq 6000 instrument (Illumina) using paired-end 150-bp reads to ~90× (n = 10) or ~30× coverage (n = 20).

To investigate how sequencing depth influenced variant detection, 10 leukemia samples (originally sequenced to 90× coverage) were down-sampled in silico to a coverage of approximately 30×. Down-sampling was performed by randomly discarding reads and retaining only every third read. The subsequent downstream analysis was performed in an identical manner to the 90× L-only analysis.

The sequence data were converted to FASTQ format using Illumina bclfastq and further processed using Sarek v2.5.1 (24), a Nextflow pipeline (25) from the nf-core framework of community-curated bioinformatics pipelines (26) on the UPPMAX Cluster (27) at Uppsala University. Briefly, preprocessed FASTQ files were checked for quality with FASTQC (28) before being aligned with BWA-mem (29) to the human reference genome build GRCh37/hg19. Duplicates were marked with Picard MarkDuplicates (30) before base calibration and indel realignment with GATK tools were performed (30), and quality statistics were aggregated with the help of multiqc (31). Following the pre-processing steps, Bam files were generated.

Sequencing data for the validation set were processed at Clinical Genomics through BALSAMIC (Bioinformatic Analysis pipeline for Somatic Mutations in Cancer) version 10.0.5, which packages the workflows for variant calling indicated below (32). Subsequently, the variant files were uploaded to the VCF visualization interface, SCOUT (33), available at https://github.com/Clinical-Genomics/scout, for further inspection and interpretation.

The vcf2cytosure v0.7.1, included in the BALSAMIC v10.0.5 converter (34) was applied to visualize ploidy changes and CNAs (35). The tool converts the output from the variant calling to a “.cgh” format used by CytoSure™ Interpret Software (Oxford Gene Technologies, Oxford, UK) originally developed to display oligonucleotide microarray measurements. Briefly, output files were binned, and 20 bins were pooled into one probe. Coverage along the genome was calculated as the log2 ratio for individual bins relative to all bins. The mean coverage in the sample was drawn at height 0 and regarded as the log2 ratio for the diploid genome (n = 2); ratios above this threshold indicate gains, and below losses. As the tool calculates the threshold relative to the average coverage over the entire genome, this threshold was adjusted manually for samples with multiple tri/tetrasomies or monosomies.

We also used allele-specific copy-number analysis of the leukemia samples (ASCAT version 4.5.0) (36), also included in BALSAMIC, to visualize aneuploidies and copy-number neutral loss of heterozygosity (CNN-LOH), with the latter requiring a paired germline sample.

Additionally, two recurrent SNVs, PAX5 P80R and IKZF1 N159Y (9), were also assessed in samples analyzed through the clinical platform.

SV calling was conducted using the FindSV (37), which merges the variants callers CNVnator (38) and TIDDIT (39), and the variant effect predictor (VEP) (40) was subsequently applied. BALSAMIC version 10.0.5 uses Manta version 1.6.0 (41), Delly version 1.0.3 (42), and TIDDIT version 3.0.0 (31) to call SVs. Subsequently, the output files were annotated in SWEGEN, a reference cohort, that reflects the genetic structure of the Swedish population (43), and variants with an allele frequency above the threshold of 0.02 were discarded. In addition to SWEGEN, the samples in the validation set were also annotated in locusDB, a locally curated database with non-tumor samples, and variants with an observed frequency above 0.02 were discarded as well (44). The SV workflow detected a median of 21,500 SVs per sample, which decreased to 7,500 once recurrent variants above the 0.02 threshold were removed (and further down to roughly 4,000 when variants in the paired germline were subtracted) (leukemia/normal, L/N). Finally, a short list of clinically relevant genes ( Supplementary Table 2 ) was applied, and the filtered events were inspected manually in IGV. The L-only analysis was carried out in the same manner, omitting the filtration with the paired normal.

The location of DUX4, within the highly repetitive D4Z4 region (45), results in a high number of copies for DUX4 and DUX4-like genes. This poses challenges for the detection of DUX4-r with short-read WGS. To overcome the limitations of our initial bioinformatics pipeline, we applied a SAMtools command (46, 47) to specifically identify reads supporting an IGH::DUX4 rearrangement. The command identifies discordant reads in the IGH region (14:106032614–107288051, in GRCh37/hg19), and subsequently, these reads were filtered to only include reads that either i) have a mate mapped to one of the regions where DUX4 may map (4:190988100–191007000, 10:135477000–135500000, or GL000228.1:7000–115000 in GRCh37/hg19) or ii) have a secondary alignment in either of the above-mentioned regions. See Supplementary Material for the SAMtools command.

The aberrations mandatory to investigate according to the current ALL treatment trial protocol (48) were all selected, as were the genes/genomic regions included in the UKALL-CNA classifier and IKZF1+ profiles (49, 50). Furthermore, the key genes signaling genetic subtypes with potential diagnostic or therapeutic implications among B-other ALL (10–16, 51–55) or included in WHO 2022 classification (8) or the ICC of hematological malignancies were also included (9). The genes’ IDs and genomic coordinates are listed in Supplementary Table 2 .

To investigate whether the detection of all mandatory genetic aberrations was feasible, 38 samples from patients with pediatric B-cell ALL, representative of the genetic subgroups stipulated by the ALLTogether trial protocol (48), were analyzed with the corresponding paired sample ( Table 1 ). These included samples with high hyperdiploidy (HeH) and mandatory aberrations caused by SVs, resulting in recurrent fusion genes or iAMP (21). Two different visualization approaches were used to visualize aneuploidies and large CNAs, vcf2cytosure and ASCAT. Both applications allowed for an accurate calling of trisomies/tetrasomies, and the findings were identical to those obtained by SoC for all 11 samples with HeH ( Supplementary Figure 1A ). Equally, the high-risk aberration iAMP (21) was readily identified by both, although inspection to determine the boundaries of the changes along chromosome 21q was only possible in vcf2cytosure ( Supplementary Figure 1B ).

Filtering SVs with the short list of genes recurrently rearranged in ALL returned 20–25 variants per sample and made manual inspection in IGV feasible. With the use of this approach, all seven samples harboring an ETV6::RUNX1 fusion were identified, irrespective of whether the fusion arose through a balanced translocation or a complex genomic rearrangement involving multiple chromosomes ( Figure 1A ). Similarly, the TCF3::PBX1 fusion, resulting from a translocation between chromosome arms 1q and 19p ( Figure 1B ), was identified in all three cases. In the initial analysis, only four out of six KMT2A-r with the corresponding fusion partner were detected ( Figure 1C ). Revision of the VCF files revealed that the missing two rearrangements had been discarded in the paired analysis with the germline. Those two follow-up samples had been taken shortly before the patients relapsed, and as the KMT2A-r already had reappeared, it had been discarded. t(9;22) was successfully detected in the two Philadelphia-positive ALL samples ( Figure 1D ), both displaying the minor BCR breakpoint and juxtaposition of exon 1 in BCR to exon 2 in ABL1. Two ABL-class rearrangements involving PDGFRB were also identified, whereas a known RANBP2::ABL1 fusion was missing from the variant list for P047. Manual scrutiny of the variants annotated by VEP in IGV revealed that the ABL1 breakpoint mapped immediately upstream of the gene. By replacing the gene’s ID with the genomic coordinates and adding 5 kb upstream, the variant was retained in the filtration step, and the RANBP2::ABL1 fusion was detected ( Figure 1E ).

Next, we assessed whether the mandatory aberrations could be confidently detected through analysis of L-only. Also in this case, filtering with the shortlist narrowed the SVs to manageable numbers and detected all the sought rearrangements, including the KMT2A-r that were not detected in the initial L/N analysis. Otherwise, the findings with both approaches were identical, and all the class-defining aberrations were detected. Supplementary Table 3 lists the findings in the samples with mandatory aberrations by SoC and WGS.

To validate the performance in the clinical setting, 10 additional samples (marked with a single asterisk, Supplementary Table 3 and Supplementary Table 4 ) were processed and analyzed through the clinical platform in a blinded manner. Subsequent comparison with SoC revealed that all relevant variants had been detected by the pipeline through the analysis of paired L/N as well as L-only. Mandatory aberrations were identified in seven samples ( Supplementary Table 3 ), two of the remaining samples harbored rearrangements of PAX5 (P099 and P100), and no single lesion was identified in the third sample (P101). The comparison to SoC showed complete concordance for all 10 cases ( Supplementary Table 4 ). In summary, WGS was able to detect all class-defining aberrations that are mandatory to investigate in the current ALL treatment trial with accuracy comparable to SoC. Figure 2A illustrates the comparison of the genetic findings by SoC and WGS for the entire set of samples.

In addition to SVs and aneuploidies, we investigated whether WGS is suitable to detect CNAs affecting the eight loci included in the UKALL-CNA classifier (49) as well as the IKZF1 ^PLUS profile (56). Deletions of various sizes from single exons to the entire IKZF1 were detected in vcf2cytosure and also by the SV callers ( Supplementary Figure 2A ), both among the established mandatory subgroups and B-other ALL. CNAs affecting PAX5, such as amplification of exons 1–5 ( Supplementary Figure 2B ), intragenic deletions, and partial deletions leading to a PAX5::ZCCHC7 fusion, were also detected in both groups, albeit at a higher frequency in B-other ALL. The same was true for losses of CDKN2A and CDKN2B. Deletions in RB1, EBF1, ETV6, ERG ( Supplementary Figure 2C ), and BTG1 were identified in isolated cases, whereas CNAs in the PAR1 region ( Supplementary Figure 2D ), which resulted in the juxtaposition of P2RY8 and CRLF2, were only observed among B-other ALL and DS-ALL. Deletion of ERG was the only lesion detected in two of the B-other samples (P050 and P114). Investigation of CNAs was not mandatory in the NOPHO 2008 trial protocol; however, the few samples with available information were concordant with the WGS findings.

A total of 36 B-other ALL and three DS-ALL, lacking stratifying aberrations in SoC, were investigated to determine the diagnostic yield. The analysis was restricted to recurrent aberrations that either define subgroups within the B-other group or have been suggested to impact outcome. Among these patients, 35 could be tentatively assigned to one of the emerging genetic categories based on the WGS findings. The lesions identified in B-other ALL are summarized in Figure 2B and Supplementary Table 4 .

The most frequent aberrations were rearrangements that affected PAX5 (n = 17) and CRLF2 (n = 10). The alterations in PAX5 (PAX5-alt) were heterogeneous and included recurrent fusions with NOL4L (n = 4) ( Supplementary Figure 3A ), ZCCHC7 (n = 3), SNTA1 (n = 1), and DACH2 (n = 1), deletions (n = 9), or amplifications (n = 2) (13). PAX5-alt often co-occurred with other putative class-defining lesions, in particular samples harboring CRLF2 rearrangements (CRLF2-r). Both previously described types of CRLF2-r were found, i.e., translocations t(X;14)(p22.33;q32) ( Supplementary Figure 3B ) juxtaposing the IGH locus to CRLF2 (n = 3) and a P2RY8::CRLF2 fusion (n = 7) (57). Interestingly, CRLF2-r always co-occurred with PAX5-r, with the exception of the three DS-ALL that carried P2RY8::CRLF2 fusions as the only detected lesion. WGS also identified four cases harboring ZNF384 fusions; the partners were TCF3 (n = 3) ( Supplementary Figure 3C ) and TAF15 in P067. An ETV6::IKZF1 fusion signaling ETV6::RUNX1-like subgroup was detected in P105. Another sample harbored an ETV6 deletion together with the KMT2C::IKZF1 fusion, and as both lesions associate with ETV6::RUNX1-like ALL, the sample was classified as such (51). A FUS::ERG fusion was detected in sample P089. This fusion is rare but recurrent in ALL and has been associated with the ETV6::RUNX1-like gene expression profile (13). Two recurrent aberrations were also found in one patient each, a MEF2D::BCL9 fusion (P078) and an IGH::ID4 fusion (P091).

In addition, WGS detected isolated ERG deletions in one sample (P050) and a rearrangement affecting the IGH locus in sample P077. The corresponding discrepant pairs mapped to the long non-coding RNA CCDC26 at 8q24; however, no cluster of discordant reads supporting DUX4-r was found in our initial analysis, despite extensive manual scrutiny in IGV. In an attempt to overcome this shortcoming, we investigated all samples sequenced to a depth of 90× (n = 68), applying a command that specifically returns the number of discordant reads that link the IGH locus to any of the copies of DUX4/DUX4-like genes in hg19. This approach was successful and identified clusters of discordant reads in these two samples and sample P101 with no lesion detected previously. As many as 126, 31, and 183 read pairs supported an IGH::DUX4 rearrangement in P050, P077, and P101, respectively ( Supplementary Table 5 ), while the median number of discordant read pairs in the remaining 65 samples was 0 (range, 0–6).

No putative class-defining lesion was identified in the remaining four B-other samples; the only lesions found by WGS were loss of CDKN2A/B, uniparental disomy (UPD) for chromosome 14, and ETV6 deletion ( Figure 2B ; Supplementary Table 4 ).

Finally, we investigated whether decreasing the coverage to 30× and L-only could be suitable in the diagnostic setting with two approaches. First, we in silico (58) down-sampled 10 leukemia samples (marked with two asterisks, Supplementary Table 3 and Supplementary Table 4 ), containing various aberrations (ETV6::RUNX1, TCF3::PBX1, KMT2A-r, BCR::ABL1, iAMP (21), and five B-other ALL) to 30× coverage and repeated the analysis described above. The comparison revealed that all the variants detected, i.e., aneuploidies, SVs, and focal CNAs, were also detectable in the corresponding down-sampled 30× data, and all 10 samples were allocated to the correct genetic subgroup.

Subsequently, DNA from 20 additional diagnostic BM samples, including two samples with low blast counts of 37% and 14%, respectively, was sequenced to 30× coverage and processed as described for L-only samples. A mandatory aberration or putative driver event was identified in 19/20 samples ( Table 2 ). Chromosomal gains as in HeH or losses as in HoL ( Figure 3A ), as well as SVs leading to recurrent fusions, were correctly detected. Even in the samples containing a low percentage of blasts, 30× coverage could identify the class-defining lesion, e.g., the ETV6::RUNX1 fusion in P111 ( Figure 3B ) and the KMT2A::AFF1 fusion P120 ( Figure 3C ). ERG deletion was the only somatic aberration identified in sample P114. Targeted analysis revealed 34 read pairs in support of an IGH::DUX4 rearrangement compared to a median of 0 (range 0 to 3) pairs for the rest of the cases, thus confirming the presence of a DUX4 rearrangement in sample P114. A summary of the findings by SoC and WGS is presented in Figure 3D . Taken together, the results indicate that decreasing the sequencing depth to 30× enables the identification of clinically relevant genomic lesions in all the samples investigated irrespective of blast count.