Integrative genomic analysis of childhood acute lymphoblastic leukaemia lacking a genetic biomarker in the UKALL2003 clinical trial

Incorporating genetics into risk-stratification for treatment of childhood B-progenitor acute lymphoblastic leukaemia (B-ALL) has contributed significantly to improved survival. In about 30% B-ALL (B-other-ALL) without well-established chromosomal changes, new genetic subtypes have recently emerged, yet their true prognostic relevance largely remains unclear. We integrated next generation sequencing (NGS): whole genome sequencing (WGS) (n = 157) and bespoke targeted NGS (t-NGS) (n = 175) (overlap n = 36), with existing genetic annotation in a representative cohort of 351 B-other-ALL patients from the childhood ALL trail, UKALL2003. PAX5alt was most frequently observed (n = 91), whereas PAX5 P80R mutations (n = 11) defined a distinct PAX5 subtype. DUX4-r subtype (n = 80) was defined by DUX4 rearrangements and/or ERG deletions. These patients had a low relapse rate and excellent survival. ETV6::RUNX1-like subtype (n = 21) was characterised by multiple abnormalities of ETV6 and IKZF1, with no reported relapses or deaths, indicating their excellent prognosis in this trial. An inferior outcome for patients with ABL-class fusions (n = 25) was confirmed. Integration of NGS into genomic profiling of B-other-ALL within a single childhood ALL trial, UKALL2003, has shown the added clinical value of NGS-based approaches, through improved accuracy in detection and classification into the range of risk stratifying genetic subtypes, while validating their prognostic significance.

We have demonstrated that whole genome sequencing (WGS) provides an excellent method for classifying B-ALL patients into clinically relevant genetic subtypes [22]. Here, we combine results from cytogenetics, fluorescence in situ hybridisation (FISH) and Multiplex Ligation-dependent Probe Amplification (MLPA) with both WGS and targeted next generation sequencing (t-NGS) of a large cohort of B-other-ALL from a single highly successful UK childhood ALL clinical trial, UKALL2003. Using this integrated approach, we have accurately classified these patients into 15 distinct genetic subtypes, described the spectrum of underlying abnormalities, and clarified their frequency and clinical significance.

Patient cohort
Patients were diagnosed with B-ALL and treated on the UKALL2003 trial (NCT00222612) (age 1-24 years) [23,24]. The Scottish Multi-Centre Research Ethics Committee approved the trial and written informed consent was obtained in accordance with the Declaration of Helsinki.
Among the total trial recruitment (n = 3204), 741 patients were classified as B-other-ALL, which excluded Down Syndrome individuals (n = 65), patients not fully tested (n = 38) and cases with normal karyotypes and ≥4 additional RUNX1 signals by FISH (n = 75), as they were likely to be undetected high hyperdiploidy or iAMP21-ALL ( Supplementary Fig. 1). Samples were available for genetic testing on a representative cohort of 351/741 B-other-ALL patients. Initially patients were assigned to three or four drug induction based on NCI risk status. High-risk patients (slow early response or high-risk genetics) were assigned to augmented postinduction therapy (regimen C), while the remaining patients were randomised to either treatment reduction (minimal residual disease (MRD) low-risk) or intensification (MRD high-risk).

Targeted NGS
Whole-genome amplification (WGA) of 30 ng genomic DNA was performed using the Repli-G Mini kit (Qiagen). RNA baits were designed to capture the whole gene sequence of 23 genes and exonic regions of 35 genes, using SureDesign and the SureSelect XT2 platform (Agilent Technologies, Santa Clara, USA), covering 97% (average) of the target regions, ranging from 69% (IKZF1) to 100% (Supplementary Table 1). DNA was fragmented into 800-1000 bp fragments by sonication using the M220 Focused Ultrasonicator (Covaris, USA) or Bioruptor Pico (Diagenode, USA). Sequencing libraries were prepared using the custom-designed SureSelect XT2 kit (Agilent Technologies, USA) according to manufacturer's protocol, with the following modifications for enrichment of larger DNA fragments: 1) 1 and 2 min annealing and elongation stage, respectively, during the pre-and post-hybridisation PCR and 2) the ratio of AMPure XP beads (Beckman Coulter, USA) to sample volume was reduced to 0.7:1. Individual samples were barcoded for pooling at equal volumes prior to sequencing. The libraries were sequenced on a HiSeq2500 or NextSeq (Illumina, UK) using 125-150 bp paired-end chemistry. Samples were sequenced to a mean coverage of 300-fold.
Raw Fastq reads were processed using the Genome Analysis Toolkit (GATK) [25]. Reads were aligned to hg19/GRCh37 and duplicates removed using BWA-MEM [26] and Picard [27]. Structural variants (SVs) were manually interrogated from deduplicated bam files, using Integrated Genomics Viewer (IGV) (Broad Institute, USA) [28], with the minimum and maximum insert size set to 50 bp and 5000 bp, respectively. Single nucleotide variations (SNVs) and indels were called using GATK HaplotypeCaller (version 3.8) [29]. Default settings were used but ploidy was increased to eight for the detection of subclonal variants and Base Quality Score Recalibration (BQSR) was not applied due to the small size of the targeted region. Hard filtering was performed and variants that passed filtering with an allele depth (AD) of 10 were annotated using Variant Effect Predictor (VEP) (version 102.0) [30], adding information from SIFT (version 5.2.2) [31] and Polyphen (version 2.2.2) [32].

Whole genome sequencing
WGS was performed on matched diagnostic and remission DNA samples, as previously described [22].

Statistical analyses
Event-free survival (EFS) was defined as time-to-relapse, second tumour or death, censoring at date of last contact. Relapse rate (RR) was defined as time-to-relapse for those achieving complete remission, censoring at date of death in remission or last contact. Overall survival (OS) was defined as time-to-death, censoring at date of last contact. The median follow-up time for the whole cohort was 10.98 years (IQR 3.83 years). Kaplan-Meier methods were used to estimate survival rates and univariate Cox regression models were used to determine hazard ratios. Other comparisons were performed using χ2 or Fisher's exact tests, as appropriate. All p-values were two-sided and values <0.05 were considered statistically significant. All analyses were performed using Intercooled Stata (Stata Statistical Software Release 16; StataCorp, USA).

Comparison of techniques
There was high concordance between WGS and t-NGS results, with the same subtype-defining abnormality identified in 28/32 (88%) cases (Supplementary Table 8). Four cases tested by both Standard-of-care techniques include cytogenetics, FISH, MLPA and SNP array. NGS includes WGS and t-NGS. Abnormal FISH signal patterns classed as balanced rearrangements: 1R1G1F, or unbalanced: 1R0G1F or 0R1G1F, with evidence of fusion from karyotype, partner gene FISH, SNP array or RT-PCR, as previously published [12]. MEF2D::CSF1R and ETV6::ABL1 are classified as ABL-class fusions, PAX5::JAK2 and ETV6::JAK2 are classified as JAK2-r, PAX5::ETV6 are classified as PAX5alt according to previously published data [2]. All other rearrangements of ETV6 are assigned to the ETV6::RUNX1-like subtype. PAX5 mutations and CN abnormalities of PAX5, CDKN2A/B and MTAP are classified as PAX5alt, only in the absence of other subtype defining abnormalities. *DUX4 and MEF2D were not included in the t-NGS kit. CN copy number, SV structural variant, CNV copy number variant.
techniques remained unclassified, although WGS identified novel abnormalities, including three in the MAP kinase pathway, as previously reported [22]. Although the t-NGS panel did not include DUX4, t-NGS identified ERG abnormalities in 8/11 DUX4-r cases analysed by WGS and t-NGS, including deletions (n = 7) and inversions (n = 1). Both WGS and t-NGS detected subtype-defining mutations in PAX5 (P80R and others) (n = 22), IKZF1 (N159Y) (n = 4), and ZEB2 (H1038R) (n = 4). They also identified secondary mutations in a wide range of genes, although none was associated with a particular subtype and numbers were small (Supplementary Table 9).
WGS detected a higher number of focal copy number abnormalities (CNA) than MLPA, notably of ERG deletions, as we have previously reported [22]. Similarly, t-NGS identified ERG deletions in three patients that had been called normal by MLPA. These deletions were either focal, with evidence of a single probe deletion only by MLPA, which is insufficient to call a deletion (n = 2), or sub-clonal, where the MLPA ratio for the deleted probes was 0.76-0.9 but not below the required 0.75 cut-off level to call a deletion (n = 1) [35,38].
FISH and MLPA detected the highest number of CRLF2 rearrangements, in 100% agreement with WGS. Detection of CRLF2-r by t-NGS was inconsistent, identifying only 2/4 IGH::CRLF2 and 2/14 P2RY8::CRLF2 cases observed by FISH and/or MLPA. This was due to incomplete coverage across the PAR1 region, with only CRLF2 being included in the t-NGS panel design, with 82% coverage.
There was high concordance between our existing FISH data and both WGS and t-NGS for other genes (98-100%, Supplementary Table 10) [12,33], with only four additional cases identified by NGS, where FISH had shown a normal result. These were MEF2D::CSF1R, SMARCA2::ZNF384, CUX1::NUTM1 and IGH::CEBPA. It is likely that these fusions resulted from complex rearrangements, as seen for MEF2D::CSF1R (9850), with multiple rearrangements, including deletions and inversions involving the MEF2D gene identified by WGS. Alternatively, the breakpoints were outside the region covered by the FISH probes, as demonstrated in the patient with the CUX1::NUTM1 fusion (20750). The translocation, t(7;15) (q22;q14), was present in the karyotype but NUTM1 was not seen to be rearranged by FISH.
FISH alone could not fully characterise rearrangements of PAX5 and ETV6. Among PAX5alt and ETV6::RUNX1-like patients, different FISH signal patterns were observed, including whole and partial deletions targeting PAX5 and ETV6, as well as balanced rearrangements (Supplementary Tables 2 & 3). Many of these cases were further characterised by t-NGS, although two cases of PAX5alt, with dic(9;20) by cytogenetics (10061 and 22861), were not called by t-NGS. Despite the presence of large abnormalities involving 9p, the breakpoints were outside the regions of PAX5, MTAP and CDKN2A/B covered by the t-NGS panel and were therefore undetected.  [36] and IKZF1 plus [40]). Coexistence of CRLF2-r is indicated in red in the B-other-ALL subtype row. Copy number profile status was unavailable for patients lacking Multiplex Ligationdependent Probe Amplification (MLPA) data. *The SALSA P335-ALL-IKZF1 and P327-iAMP21-ERG MLPA kits were used to determine gene copy number. Relapses were defined as follows: very early, < 18 months from diagnosis; early, within 6 months of end of treatment; late >6 months after end of treatment. B Circos plot illustrating the range of PAX5 translocation partner genes in PAX5-r. C Stacked bar plot showing the distribution of PAX5alt abnormalities amongst different age groups. NCI National Cancer Institute, WBC White blood cell count, MRD Minimal residual disease, t-NGS Targeted next-generation sequencing, WGS Whole genome sequencing, PAX5-r PAX5 rearranged, PAX5-ITD PAX5 internal tandem duplication, CNA Copy number alteration.

B A
while ERG deletions were identified in a further 19 cases tested by t-NGS and/or MLPA ( Supplementary Fig. 4).
Rearrangements of CRLF2 were observed in 53 patients, with P2RY8 (n = 33) or IGH (n = 20) partners, including the nine cases mentioned above, where the fusion co-existed alongside other subtype-defining abnormalities.
ZNF384 fusions were found in 37 patients, involving nine different partner genes, including SPI1, which has not previously been reported in B-ALL ( Supplementary Fig. 6).
New genomic subtypes and additional risk factors CNA identified from the P335 MLPA kit varied according to subtype, with several associations reaching statistical significance (Table 3). Previous studies have reported that specific copy number profiles have prognostic relevance. We investigated the distribution of the copy number profiles, UKALL-CNA [36,39] or IKZF1 plus [40] . among B-other-ALL (Table 3). Patients with PAX5alt and CRLF2-r had an increased frequency of the poor-risk IKZF1 plus profile and were more likely to be classified as UKALL-CNA intermediate/poor-risk (IR/PR) compared to patients in other subtypes. ZNF384-r and DUX4-r cases were more likely to have a UKALL-CNA good risk (GR) profile. Unsurprisingly, patients with DUX4-r were less likely to have an IKZF1 plus profile, given the association with ERG deletions.
Where numbers permitted, we assessed whether the presence of deletions or copy number profiles modulated the outcome of patients within subtypes (Supplementary Table 11). Survival rates for patients with PAX5alt and CRLF2-r, who also had an IKZF1 deletion or IKZF1 plus profile, appeared inferior, although log rank tests revealed borderline p values, suggesting that they were not the main drivers of poor outcome. In contrast, at end of induction (EOI), high levels of MRD were strongly associated with increased RR and lower EFS within the PAX5alt subtype. The UKALL-CNA profile was too tightly correlated with many subtypes to be assessable but was linked to outcome in ZNF384-r cases. Further analysis of PAX5alt and CRLF2-r revealed no difference in outcome within each subtype, according to NCI risk group.

DISCUSSION
In this study, we have comprehensively refined the classification of B-other-ALL by integrating NGS-based techniques with those that we have previously reported [12,33]. We have demonstrated the value of incorporating both WGS and t-NGS for improved identification of a range of abnormalities associated with emerging subtypes, particularly, for detection of subtypedefining mutations.
While our previous approaches were highly successful in classification of B-other-ALL, sequencing-based methods provided valuable additional information in many cases. The increased sensitivity of NGS identified the full range of secondary and cooperating abnormalities, for example, ERG abnormalities in DUX4-r patients. Notably, NGS identified fusion partners, whereas FISH detected only the rearrangement of the relevant "hub" gene, such as ZNF384 or PDGFRB. These data may be important in future collaborative studies, to discern clinical associations for specific fusion genes within subtypes. For example, a recent international collaboration collected data from 218 patients with ZNF384-r and showed EP300::ZNF384 to be associated with a lower risk of relapse compared to other ZNF384 fusions [44].
NGS approaches were particularly informative in defining the genomic abnormalities characteristic of two subtypes, PAX5alt and ETV6::RUNX1-like, previously identified from gene expression profiling. Both subtypes are associated with a variety of abnormalities, which differ between patients and occasionally overlap with other subtypes, rendering them difficult to define by standard-of-care techniques. Building on our recent WGS study [22], here we have demonstrated that NGS can reliably detect these subtypes prospectively within a diagnostic setting without the need for expression profiling.
Neither DUX4 nor MEF2D were included in the t-NGS kit, as they were unknown at the time of design, thus highlighting the importance of flexibility when choosing tools for diagnostic testing. We were able to screen for MEF2D rearrangements by FISH [12], however, accurate DUX4-r identification was only possible using WGS with a bespoke analysis pipeline [22]. It remains to be determined whether standard PCR testing or t-NGS with a similar customised pipeline can reliably identify DUX4-r. In our recent WGS study, the occurrence of an ERG abnormality was pathognomonic of the DUX4-r subtype, although only present in 68% of cases [22], thus reliance on ERG deletion detection as a surrogate marker of the DUX4-r subtype would miss > 30% cases.
Due to the relatively small numbers of previously published cases, the true prognostic impact of these subtypes remains unresolved. The excellent long-term survival for DUX4-r patients in this study has extended the observations made by others, reporting high 5-year survival rates [3,41,42]. There is growing evidence that DUX4-r patients have low RR; possibly linked to the increased therapy that they receive based on EOI MRD [41,42]. Compared to patients with ETV6::RUNX1 and high hyperdiploidy, DUX4-r patients were more often NCI and MRD high-risk, so more likely to be treated on more intensive treatment regimens (Supplementary Table 12). This phenomenon was mentioned in previous DUX4-r/ERG deletion studies [3,41,43,45,46]  question as to whether their excellent outcome was due to intensified treatment or that DUX4-r is an intrinsically chemosensitive good-risk subtype. Here we have shown no evidence that relapse is linked to therapy. Moreover, due to their long-term excellent outcome [47,48], it is reasonable to consider these patients as cured.
It is now widely recognised that patients with ABL class-fusions not treated with tyrosine kinase inhibitors have a very poor prognosis [21,49], as further reinforced here. No relapses or deaths were observed among 21 patients with ETV6::RUNX1-like-ALL after a median follow-up of 10 years. This excellent outcome differs from 22 and 13% 5-year cumulative incidence of relapse reported for 18 ETV6::RUNX1-like-ALL cases treated on Total Therapy 16 (n = 9) [42] or MS2003/2010 (n = 9) [41], respectively. Although here we primarily used DNA-based techniques to identify the genomic abnormalities associated with this subtype, we confirmed an ETV6::RUNX1-like gene expression signature in six patients by WTS [22]. Only four of the 21 ETV6::RUNX1-like patients received intensive therapy, suggesting that, when treated on UKALL2003, they have an excellent outcome. The outcome of patients with PAX5alt, CRLF2-r or ZNF384-r was very similar to B-other-ALL overall, broadly consistent with other paediatric ALL trial publications [41,42,50]. The MS2003/2010 study reported an adverse effect of IKZF1 deletions within the PAX5alt group [41]. Although our results were consistent with their observation, it was eclipsed by the negative effect of MRD. We identified too few patients with PAX5 P80R, IGH::ID4, ZEB2/CEBP, MEF2D-r, NUTM1-r or IKZF1 N159Y to reliably assess outcome. Given their rarity, international collaborative studies are needed to determine their true risk status.
It is evident that accurate classification of B-other-ALL is crucial to the success of future trials, thus access to a range of approaches for their detection is important. Other studies have applied WTS and subsequent cluster analysis to retrospectively classify B-other-ALL [41,42,51]. Here, we have chosen DNA-based approaches, for detection of the defining genetic abnormalities. As our associated study demonstrated high concordance between WGS and WTS, we are confident that our genomic approaches, specifically WGS, can accurately and prospectively classify B-other-ALL [22].
Our methodology has a number of advantages: while WTS requires a large reference cohort, these DNA-based techniques can be performed on individual or small numbers of cases. Both WGS The UKALL-CNA or IKZF1 plus profiles are based on the genes included within the P335-IKZF1 MLPA kit. Briefly, the UKALL-CNA profile classifies patients as good risk (CNA-GR), if they have no deletions among the genes tested for, isolated deletions of ETV6, PAX5, BTG1 or ETV6 with a single additional deletion of BTG1, PAX5, or CDKN2A/B. All other profile combinations are classified as intermediate/poor-risk (CNA IR/PR) [36]. The IKZF1 plus profile defines patients with an IKZF1 deletion and at least one additional deletion of PAX5, CDKN2A/B, or PAR1, in the absence of an ERG deletion, as poor-risk [40].
and WTS are costly, requiring sophisticated bioinformatics pipelines for analysis, which will be prohibitive for many low and middleincome countries. As this study has demonstrated a high level of concordance between WGS and both t-NGS and standard techniques, although developed countries may adopt WGS as the predominant diagnostic test for ALL in future, laboratories with limited resources may choose standard techniques to screen only for those abnormalities linked to treatment implications. For example, in UK trials, we have previously shown that FISH testing for ABL-class fusions in patients with refractory ALL is highly effective for identification of the majority of patients [21,52]. Choice of diagnostic testing will also be driven by the preferences of different centres and be dependent on individual trial requirements.
In conclusion, we have successfully classified 351 patients with B-other-ALL into key genomic subtypes, using both NGS and standard techniques; thereby providing screening options to suit all resource levels and trial protocols. As this study was based on a single clinical trial, we were able to provide robust and clinically useful prognostic information on six recently reported genomic subtypes.

DATA AVAILABILITY
DNA and RNA Sequencing data have been deposited in the European Genomephenome Archive (EGA) under the Accession Code EGAS00001006863. Alternatively these data will be made available upon request from Dr. Sarra Ryan (sarra.ryan@newcastle.ac.uk) or Prof Christine Harrison (christine.harrison@newcastle.ac.uk).