Classification and genetics of pediatric B-other acute lymphoblastic leukemia by targeted RNA sequencing

Acute lymphoblastic leukemia (ALL) is the most common childhood malignancy, accounting for 25% of pediatric cancers.¹ Over the last decades, advances in ALL treatment and supportive care have dramatically increased the survival rate of pediatric ALL to between ∼80% and 90%.^2-4 Contemporary protocols rely on sophisticated risk stratification of patients, which includes genetic features^5,6 and minimal residual disease levels.^7,8 In this context, the genetic risk factors include recurrent chromosomal rearrangements, aneuploidy, and cooperating mutations involved in leukemogenesis.

Approximately 70% of B-cell precursor ALL (BCP-ALL) can be classified by standard genetic analyses into well-known subtypes, based on the presence of t(9;22) BCR-ABL1 (Ph-positive), rearrangements of KMT2A (MLL) at 11q23, t(1;19) TCF3 (E2A)-PBX1, t(17;19) TCF3-HLF, t(12;21) ETV6-RUNX1, high-hyperdiploidy (HeH), or hypodiploidy.^9,10 Until some years ago, the remaining 30% of pediatric BCP-ALL were called “B-other”, because this portion did not have any of the aforementioned characteristic chromosomal aberrations at diagnosis, and the underlying driver events were unknown.⁹ After extensive genetic studies, some B-other ALL cases could be grouped into common; nonoverlapping genetic subgroups, such as Ph-like, iAMP21, DUX4-cluster, ZNF384-rearranged, MEF2D-rearranged, and NUTM1-rearranged; and PAX5-driven.^11-13 Identifying these new biomarkers are important because some serve as novel therapeutic targets and/or provide prognostic information.¹⁴ 

Current ALL diagnostic strategies include karyotyping or DNA index by cytometry to detect high-hyperdiploidy, and reverse transcription polymerase chain reaction (RT-PCR) or fluorescence in situ hybridization (FISH) for “classical” fusions. Patients at high risk who are negative for this initial testing are eligible for RT-qPCR, FISH, Sanger and/or next generation sequencing (NGS) to detect other actionable/prognostic biomarkers, like CRLF2-high, JAK1/JAK2 mutations, and ABL-class fusions.^15,16 This complex workflow might be simplified by NGS methods and extended to all patients.

Here, we present a study on 144 pediatric B-other ALL cases in conjunction with 40 BCP-ALL cases of the classical genetic groups, to assess the feasibility and benefits of introducing a commercial panel of targeted RNA sequencing into the routine diagnostics. Moreover, we present, to the best of our knowledge, for the first time, the frequencies of these novel ALL subtypes in a South American population and explore their associations with clinico-biological features at presentation.

This study included 144 of 177 consecutive cases of B-other ALL treated in a single institution and 40 cases of classical ALL subtypes included as controls: high-hyperdiploidy (n = 3), TCF3-PBX1 (n = 3), ETV6-RUNX1 (n = 24), BCR-ABL1 (n = 4), and KMT2A-rearrangements (n = 6). Patients (≤18 years) were diagnosed and treated at the Boldrini Children’s Center in accordance with the Brazilian GBTLI ALL-1999¹⁷ and -2009¹⁸ protocols. General demographic and genetic information of patients enrolled in these protocols is summarized in supplemental Table 1. Leukemia diagnosis was established based on morphology, immunophenotype, cytogenetics, RT-PCR, and multiplex ligation probe-dependent amplification (MLPA, using kits from MRC Holland) as part of routine clinical diagnostics. B-other classification was established by a lack of t(12;21), t(1;19), t(9:22), and KMT2A-rearrangements and by having ≥44 and <51 chromosomes (or DNA index of <1.16). The study was approved by the institutional ethical committee (CAAE51991015.0.0000.5376). All patients, parents, or guardians provided informed consent for sample collection and research in accordance with the Declaration of Helsinki.

RNA was obtained from bone marrow mononuclear cells. RNA integrity (DV200 > 50%) was checked on an Agilent 2100 Bioanalyzer (Agilent Technologies). Targeted RNA sequencing of 1385 cancer-related genes was done using the TruSight RNA Pan-Cancer kit (Illumina) per the manufacturer’s instruction. Approximately 3 million reads per sample (2 × 76 cycles) were sequenced in MiSeq equipment.

Data analysis was performed using Illumina BaseSpace apps: RNA-seq Alignment version 1.1 and version 2.0, and DRAGEN RNA Pipeline version 3.5, with standard settings. Reads were mapped to the GRCh37/hg19 reference genome. Fusion detection was also performed using CICERO.¹⁹ Fusions were considered true positive if reported as “PASS” by at least 2 software. IGH, ERG, and CRLF2 fusions were considered true positive even by a single software (usually DRAGEN RNA or CICERO). Predicted fusions were manually curated on ProteinPaint/FusionEditor (https://proteinpaint.stjude.org/FusionEditor/), excluding fusions caused by misalignment of reads to pseudogenes. All rearrangements were manually reviewed using the integrative genomics viewer (IGV).²⁰ Detection of PAX5 amplifications, CRLF2-rearrangements, and ERG and IKZF1 intragenic deletions were manually curated by analyzing “soft-clipped” and exon junction reads in IGV. Strelka variant caller was used for single-nucleotide variant (SNV) and small indel identification. Salmon was used for gene-expression quantification in transcripts per million (TPM). Strelka and Salmon were run in the RNA-seq Alignment version 2.0 platform. A subset of genomic data for representative cases of each of the alterations described in our study was submitted on Zenodo under URL https://doi.org/10.5281/zenodo.7559347.

Individual gene-expression values (TPM) were log2 transformed and then subjected to Euclidean hierarchical clustering analysis (HCA) with a method of Ward. Different probe sets used for ALL clustering analysis were selected from the literature and genes matching with the 1385 targeted panel were used for HCA. The final lists of genes used in clustering analysis are provided in supplemental Table 2. For t-distributed stochastic neighbor embedding (t-SNE) analysis, the log2 TPM expression values were filtered using the varFilter function on the genefilter R package (version 1.74.0), with a cutoff set at 0.5 and with interquartile range (IQR) as the base function. The resulting 731 most variable genes were used for t-SNE. The t-SNE plot was generated using the tsne function on the M3C R package (version 1.14.0)²¹ using a perplexity score of 10 (set.seed(1238) perplex = 10).

The rs3781093 and rs3824662 GATA3 polymorphisms were genotyped by a 2-step PCR (primers in supplemental Table 3) followed by Sanger sequencing.²² For GATA3 allele frequencies in non-ALL controls, we used the ABraOM (Brazilian Genomic Variants) database, consisting of 1171 admixed individuals from São Paulo (Brazil).²³ The validation of fusion transcripts was done by RT-PCR and/or FISH, by routine laboratory procedures.

Possible associations between genetic abnormalities and demographic or clinical data were estimated using a 2-sided Fisher exact test with Monte Carlo simulations (based on 10 000 sampled tables with starting seed of 2 000 000). Analyses were performed using the SPSS Statistics, WinSTAT, GraphPad Prism, and/or R (www.r-project.org) software.

In total, ∼32% (177 of 561) of non–Down syndrome (DS) BCP-ALL cases treated at the Boldrini Children’s Center (Brazil) from 1999 to 2017 could not be classified in any of the classical molecular subtypes of ALL and were provisionally classified as B-other ALLs (supplemental Table 1). Of these, 144 cases had frozen samples available and were selected for the study. The median age at diagnosis was 7.38 years, males were predominant (55.6%), and 86% were standard risk (SR) at diagnosis based on the National Cancer Institute (NCI) risk classification (Table 1). Of the 144 cases, 7 had intrachromosomal amplification of chromosome 21 (iAMP21) by FISH and MLPA (supplemental Figure 1). We kept iAMP21 as part of B-other ALLs, because it was still a provisional entity when this work began.¹⁰ 

Except for the high-hyperdiploid cases, all other classical subtypes of ALL included as controls were correctly detected by the Pan-Cancer panel. Fusion transcripts were found in 32.6% (47 of 144) of B-other samples investigated (Figure 1). All the recently described B-other fusion transcripts were detected, including ABL-class (n = 7, 4.9%), JAK2-rearranged (n = 3, 2.1%), P2RY8-CRLF2 (n = 6, 4.2%; 1 of which was only found by manually searching soft-clipped reads in the IGV software), ZNF384-rearranged (n = 9, 6.3%), MEF2D-rearranged (n = 4, 2.8%), NUTM1-rearranged (n = 4, 2.8%), and ETV6-rearranged (n = 2, 1.4%) (supplemental Figure 2-3). Commonly missed by many popular fusion detection tools, we found 1 IGH-EPOR fusion but only by using CICERO (supplemental Table 4). IGH-CRLF2 cases could not be detected in fusion transcript analysis but were later identified by gene-expression analysis (further sections).

We found novel gene partners for JAK2 (TOP2B-JAK2), ZNF384 (NCOA3-ZNF384, SPI1-ZNF384), MEF2D (MEF2D-PYGO2), NUTM1 (KAT6A-NUTM1), ETV6 (ETV6-MANSC1), and PAX5 (PAX5-ZBED3) in the B-other cohort. Other novel in-frame fusions are shown in Figure 1C; supplemental Table 4, a subset of which were validated by FISH, RT-PCR, or MLPA (supplemental Figure 2).

There are 2 subtypes of B-ALL defined by PAX5 alterations: PAX5-altered (PAX5alt) and PAX5 p.Pro80Arg mutated (PAX5 P80R). The former comprises diverse PAX5 rearrangements, including fusions, focal intragenic amplifications, and mutations (other than P80R). Fusion transcript analysis resulted in the identification of 11 cases with PAX5 fusions, 4 out-of-frame (supplemental Table 4). Two out-of-frame PAX5-rearranged cases (PAX5-ZCCHC7, PAX5-MLLT3) also harbored a P2RY8-CRLF2 fusion. One of the ETV6-RUNX1 controls also carried a PAX5-ZCCHC7 fusion. To detect PAX5 intragenic amplification, soft-clipped sequences were manually inspected in the IGV software, as previously described.¹³ Four patients with B-other ALL showed altered PAX5 transcripts indicative of intragenic exons 2-5 amplification (supplemental Figure 4), which was confirmed by MLPA (data not shown). SNVs in the PAX5 coding and exon/intron regions were found in 12 cases, 5 of which corresponding to p.P80R. The variant allele frequency of PAX5 P80R for these patients was 44%, 49%, 72%, 99%, and 100%, respectively (messenger RNA level). Loss of the other PAX5 allele is a hallmark of PAX5 P80R.¹² The first patient harbored a frame-shift mutation in the other PAX5 allele whereas the last 2 patients had a single PAX5 loss by MLPA (data not shown). In total, we identified 13 PAX5alt (9%) and 5 PAX5 P80R (3.5%) cases.

Five B-other ALLs had a fusion between ERG and LINC01423, likely resulting from an intrachromosomal deletion at the ERG locus (ERGdel). When ERG soft-clipped sequences were manually searched in IGV, we were able to detect another 7 cases with ERG-LINC01423 fusions, totaling 12 cases. This ERG fusions seemed characteristic of some ERGdel and were easily visualized in a Sashimi plot²⁴ as splice junction reads connecting ERG exon 4 (or exon 11 in few cases) to LINC01423 (supplemental Figure 5). One of these ERG-LINC01423 cases harbored a NUP214-ABL1 fusion whereas all the other 11 were later shown to belong to DUX4-rearranged (see further sections). However, when MLPA was used to investigate ERG deletions in the entire cohort (supplemental Table 5), only 9 of 12 ERG-LINC01423 were confirmed. Moreover, MLPA detected 7 additional cases (5 belonging to DUX4-rearranged and 2 to other groups; supplemental Figure 5), not found by Pan-Cancer RNA sequencing. Therefore, two-thirds of all ERG deletions could be found by Pan-Cancer RNA sequencing. Interestingly, MLPA analysis revealed that ERG-LINC01423 was associated with ERG exons 5-11 deletions, an alteration found exclusively in DUX4-rearranged cases.

MLPA analysis revealed that 30% (41 of 136) of cases bore different IKZF1 deletions (IKZF1del). Having this information, we investigated whether IKZF1del could be identified by splice junctions track analysis of soft-clipped sequences in the IGV software. The IKZF1 gene has a relatively complex pattern of alternative splicing, making Sashimi plots hard to interpret. Nevertheless, by knowing beforehand which IKZF1 exons were deleted, we were able to find certain exon-exon connecting arcs, with at least 4 reads, unique to intragenic IKZF1del cases (Figure 2; supplemental Figure 6A). For example, the second most common type of IKZF1del, that is, exon 4-7 deletions (Δ4-7), were characterized by >200 reads spanning exon 3-8 junctions (Figure 2B; supplemental Figure 6B). The same 3-8 junctions, but with less reads, were seen in a case with Δ5-6. IKZF1 Δ4-8 showed junction reads spanning IKZF1 exon 3 to downstream genes, FIGNL1 or DDC (Figure 2C; supplemental Figure 6C).

We noticed that junctions between exons 1-5 or 1-8 and junctions between IKZF1 exon 3 and FIGNL1 did not occur in nondeleted samples, being exclusive to IKZF1del cases. Therefore, we believe these 3 junctions can be trusted to discriminate some IKZF1del cases. In summary, only the following IKZF1del isoforms could be realistically detected: Δ1-3, Δ1-7, Δ2-3, Δ2-7, Δ4-7, Δ4-8, and Δ5-6, which represent approximately half of IKZF1del cases. Neither whole gene IKZF1 deletions (Δ1-8), which account for ∼40% of IKZF1del cases, nor isoforms Δ1 and Δ1-6 could be detected by splice junction read analysis (supplemental Table 6). Representative Sashimi plots for all different IKZF1del are presented in supplemental Figure 6.

In addition to identifying recurrent fusion genes, targeted RNA sequencing simultaneously measures the expression of all captured genes within each sample. Ranking of CRLF2 expression values (log2 TPM+1) in all samples resulted in a sigmoid curve showing a group of 11 top outliers with CRLF2 levels 3 times higher than the median (Figure 3A). As expected, all P2RY8-CRLF2 cases found by fusion transcript analysis were in the CRLF2-high group (supplemental Figure 7A-7B). Four of the 5 remaining outliers were found to bear an IGH-CRLF2 rearrangement by FISH (supplemental Figure 7C). Unfortunately, there was no material for FISH analysis of the last case, which, however, had CRLF2 p.F232C mutation commonly seen in CRLF2-high cases.²⁵ The frequency of CRLF2 rearrangements was consistent with the 1:1 proportion of translocations (IGH-CRLF2) to deletions (P2RY8-CRLF2) usually found in non-DS pediatric ALL.²⁶ Of note, patients in the CRLF2-high group also had JAK2 (p.R683, n = 3; p.A658S, n = 1; p.L925P, n = 1) or CRLF2 (p.F232C, n = 3) mutations that are known to be associated with induction of ligand-independent activation of downstream signaling pathways.^27,28 A similar approach was used to analyze EPOR expression. A single case of EPOR-high was found (Figure 3B), which corresponded to the IGH-EPOR detected by fusion transcript analysis.

DUX4-rearranged ALL originates from translocation of a truncated copy of DUX4 into the IGH or ERG (less common) locus, resulting in a chimeric transcript formed by most of the DUX4 coding sequence plus noncoding regions of IGH or ERG.²⁹ We could not detect DUX4-rearrangements by fusion transcript analysis. However, by simply analyzing read coverage over DUX4 loci in chromosomes 4 and 10, and the ERG locus in chromosome 21, we differentiated a total of 35 cases showing high DUX4 gene family expression, of which 27 expressed also an alternative isoform of ERG (ERGalt) with an extra exon 6 (Figure 4). As expected, 11 of 12 ERG-LINC01423 cases (see earlier results) were among the 35 cases in the DUX4-rearranged subgroup.³⁰ 

Because the subtypes of ALL are characterized by distinct gene-expression profiles, we investigated whether unsupervised clustering analysis would differentiate the subtypes of BCP-ALL. Pan-Cancer genes matching genes previously used for ALL clustering (supplemental Table 2) were used for HCA. As shown in Figure 5A, using 124 genes of the 864 selected by Liu et al (2016),¹¹ there was a clustering tendency for some ALL subtypes, particularly DUX4-, ZNF384-, and MEF2D-rearranged ALL. Some groups were better clustered when using a different gene set, although segregation was never perfect (supplemental Figure 8). Curiously, HCA based on 45 genes differentially expressed by Ph-like ALL (supplemental Table 2) resulted in a group of 9 of 11 ALL cases carrying ABL or JAK/STAT-activating fusions and 5 of 11 CRLF2-rearranged cases (2 with P2RY8-CRLF2, 2 with IGH-CRLF2, and 1 with CRLF2 F232C mutation) (supplemental Figure 8A). The fact that half of CRLF2-high and almost all ABL and JAK/STAT-activating rearrangements were in the same group was consistent with Ph-like ALL.¹⁴ 

Importantly, when gene-expression profiles were analyzed by t-SNE, DUX4-rearranged and ETV6-RUNX1 cases were clearly segregated (Figure 5B). Besides, ZNF384-rearranged, ABL or JAK/STAT-activating fusions, iAMP21, and PAX5-driven cases were somehow grouped but mostly intermingled, which may be a consequence of the low number of cases included in the analysis.

In summary, after identifying all the group-defining genetic alterations, a total of 96 of 144 B-other cases could be further classified into 1 of the modern subgroups of ALL. For the remaining 48 cases (33.3%) (B-“rest”), no new recurrent biomarker could be found.

After demonstrating that the Pan-Cancer assay was able to identify several of the new subgroups of B-other ALL, we searched for possible associations between these subgroups and some of the patients’ features at presentation (Figure 6A). After testing for gender biases, we observed that girls were significantly enriched in the B-rest subgroup (P = .049), that is, 48 B-other cases without any known genetic marker. PAX5-driven cases tended to occur more frequently in boys (Fisher exact test; P = .043; Monte Carlo; P = .071) (Figure 6B).

To increase the statistical power, ABL-, JAK2-, and EPOR-fusion ALLs were analyzed as a single group named AJ2E-fusion. Both AJ2E (P = .007) and CRLF2-high (P = .042), which are characterized by an activated kinase signature, were significantly associated with higher white blood cell (WBC) counts at diagnosis. Moreover, the AJ2E-fusion group was associated with NCI high-risk classification (P = .038) and with higher frequency of IKZF1del (P < .0001) (Figure 6D; supplemental Table 7).

In total, 8 of the 144 B-other ALLs had CALLA (CD10)-negative immunophenotype, which were significantly enriched in the ZNF384-rearranged group (P < 10^-7). Noteworthy, 2 of 4 NUTM1-rearranged cases were infants, confirming the notion that NUTM1-rearrangements are rare but recurrent in infants, particularly among those who lack KMT2A-rearrangements (Figure 6D; supplemental Table 7).^31,32 No other associations with clinical and demographic data were observed (supplemental Table 7).

Independent genome-wide association studies in different ethnic populations have identified several risk loci for ALL susceptibility, including ARID5B, IKZF1, CEBPE, PIP4K2A, CDKN2A, and GATA3. Some of these polymorphism were associated with distinct ALL subtypes.^22,33 Notably, 2 inherited genetic variations in GATA3, single-nucleotide polymorphisms rs3824662 and rs3781093, are strongly associated with leukemia treatment outcomes and are particularly relevant in patients of Hispanic ancestry.³⁴ The susceptibility related to GATA3 polymorphisms appears to be affected by cytokine and kinase-activating alterations, as in cases of Ph-like ALL.³⁵ We found the A risk allele at rs3824662 and the C risk allele at rs3781093 to be overrepresented in most B-other ALL subgroups, particularly in CRLF2-high (Figure 6C). Interestingly, B-rest had GATA3 allele frequencies similar to that of healthy individuals.^22,23 

The laboratory methods used in the clinical diagnosis of genetic alterations in BCP-ALL include karyotyping, FISH, and RT-PCR. These methods restrict the number of genetic loci that can simultaneously be analyzed and are limited to the classical subtypes of ALL. NGS RNA sequencing circumvents this problem by providing the analysis of several loci at once and is particularly useful because it allows the identification of fusion transcripts derived from recurrent chromosomal translocations that are the hallmarks of ALL.^36-46 However, in the current diagnostics of ALL, NGS is used in addition to previous methods in intricate workflows. Our results suggest that targeted RNA sequencing might substitute FISH and RT-PCR methods, simplifying the current diagnostic strategy for ALL, while providing identification of both classical and modern subgroups of ALL, except for ploidy and copy number alterations. Despite the relatively low number of positive controls analyzed, fusion detection by targeted RNA sequencing was highly consistent with results obtained by cytogenetic and RT-PCR for BCR-ABL1, ETV6-RUNX1, TCF3-PBX1, KMT2A-AFF1, KMT2A-MLLT1, and KMT2A-MLLT3. Except for iAMP21 cases, targeted RNA sequencing allowed detection of all other novel subgroups of B-other ALL, which includes CRLF2-high, ABL-fusion, JAK2-fusion, EPOR-fusion, ZNF384-rearranged, NUTM1-rearranged, MEF2D-rearranged, ETV6-fusion, PAX5-altered, PAX5 P80R, and DUX4-rearranged (including segregation of cases with ERGalt).

Targeted RNA sequencing takes longer than total RNA sequencing but is less expensive, provides a deeper coverage of transcripts, and is easier to analyze. Moreover, targeted RNA sequencing is applicable to lower quantities/quality of RNA. We successfully used the TruSight Pan-Cancer kit for the analysis of bone marrow smears (data not shown). The advantage of using a commercial sequencing assay, with global kit distribution, is that analytical performance has been evaluated by the manufacturer and by independent researchers. The TruSight series of sequencing panels ensure deep, uniform coverage, which is 1 of the most important parameters for highly accurate variant calling. In hematologic malignancies, the TruSight Pan-Cancer was shown to detect 100% of fusions transcripts.^42,43 Regarding SNV, concordance between Pan-Cancer RNA sequencing and exome analysis was 86%, except for nonsense mutations likely subjected to messenger RNA decay.⁴² Serial dilution of BCR-ABL1⁺ RNA into normal control RNA revealed a detection limit ranging from 1% to 10%.⁴⁷ 

The TruSight RNA Pan-Cancer kit has been tested with success within the diagnostic workflow of the international AIEOP-BFM ALL 2017 trial in a cohort of 117 BCP-ALL cases.⁴³ All classical group-defining fusion transcripts were identified as well as some of the rearrangements involving CRLF2, PAX5, JAK2, PDGFRB, ETV6, ZNF384, MEF2D, and NUTM1 genes, commonly found in B-other ALL. However, the number of B-other cases included in this previous study (n = 36) was relatively small. In this manuscript, we used the Pan-Cancer panel to analyze 144 B-other ALL in conjunction with 40 BCP-ALL of the classical groups as controls. Besides confirming previous observations on the utility of this targeted RNA-sequencing panel in identifying group-defining fusion transcript, we also explored SNVs, exon junction tracks, and gene-expression data (see workflow in supplemental Figure 9), to show for the first time its utility in the identification of some other molecular markers of ALL, particularly PAX5 P80R, IGH-CRLF2, intragenic ERGdel and IKZF1del, and DUX4-rearrangements. The frequencies of these genetic alterations in the Brazilian cohort of B-other ALLs (Figure 7) seemed compatible to those found by total RNA sequencing in a recent European study of unselected consecutive B-other cases,¹³ except that we found a much higher number of fusions resulting in ABL and JAK/STAT kinase activation (11 of 144 compared with 1 of 110). We believe that this discrepancy is not a consequence of sampling bias because our inclusion rate was of 86.1% (supplemental Table 1). The discrepancy might relate to methodological differences: targeted RNA sequencing vs total RNA sequencing, read depth, and/or bioinformatics tools. Another possible reason for the increased fraction of ABL- and JAK-STAT-activating fusions in the South American cohort could be related to ancestry.⁴⁸ 

A series of novel fusion partner genes were identified: TOP2B for JAK2, KAT6A and ZNF618 for NUTM1, MANSC1 for ETV6, NCOA3 and SPI1 for ZNF384, and PYGO2 for MEF2D. Although the Pan-Cancer targets only 1385 genes, new in-frame fusion genes (JSRP1; OAZ1-DOT1L, UBE2D3-AFF1, PPFIBP1-KDM5A, and ERBB4-IKZF2) were identified in ALL cases with no other known driver alteration, deserving further studies. All these genes were already described as participating in fusion transcripts in the FusionGDB2 database (https://compbio.uth.edu/FusionGDB2/). However, to the best of our knowledge, the combinations found here have not yet been reported, except for IKZF2-ERBB4 fusion that was found in T-cell lymphomas⁴⁹ and ovarian tumors.⁵⁰ 

Importantly, some of the subgroup-defining fusion genes could not be found by fusion transcript analysis. Fortunately, the identification of these cases was still possible by means of analyzing gene-expression levels and read coverage. Patients with IGH-CRLF2 or IGH-EPOR translocations express exceptionally high levels of CRLF2 or EPOR, respectively. NUTM1-fusion cases also had high levels of expression of NUTM1, a gene that is normally not expressed in ALL (supplemental Figure 3F). Likewise, in comparison with other ALLs, DUX4-rearranged ALL showed a significantly higher pattern of read coverage along the DUX4 gene cluster. The DUX4 subgroup is enriched in ERG deregulation, that is, aberrant ERG expression, expression of an alternative splicing isoform of ERG (ERGalt), or ERG deletions.³⁰ ERGalt was identified in 27 of 35 cases of DUX4-rearranged ALL. Finally, DUX4-rearranged ALL formed an independent group in hierarchical clustering analysis, and both DUX4-rearranged and ETV6-RUNX1 ALL were easily segregated in t-SNE analysis. HCA was tentatively used to identify Ph-like ALL but with ambiguous results. Because targeted RNA sequencing identified all the ABL-class and JAK/STAT-activating gene alterations, which is decisive in guiding the choice of an alternative drug use, there is probably no need for Ph-like classification.

Focal IKZF1del are frequently found in B-other ALL. Although IKZF1del does not define a separated subgroup of ALL, identification of IKZF1del is of high relevance as a prognostic factor, being associated with worse outcome.^51,52 Here, we showed that certain exon-exon junctions captured by targeted RNA sequencing allowed the discrimination of some intragenic IKZF1del using the Sashimi plot tool in the IGV browser. However, this method is not straightforward. Moreover, whole IKZF1 gene deletions (exons 1-8), which account for one-third of IKZF1del cases could not be discriminated. Thus, Pan-Cancer RNA sequencing cannot be a substitute for current MLPA or microarray methods, which are necessary to analyze IKZF1del and other copy number alterations of prognostic significance.

Although B-other ALL is more frequent in males,¹³ we found, to the best of our knowledge, for the first time, that PAX5-driven B-other ALL and B-rest ALL groups had a significantly higher proportion of males and females, respectively. We have no explanation for these gender biases. As expected, ZNF384-rearrangements were associated with lack of CD10⁵³ and CRLF2-high with an initial WBC count of ≥50 × 10³/μL and GATA3 risk alleles (rs3781093 and rs3824662).²² The group of ABL/JAK2/EPOR-fusions was associated with high initial WBC counts, high risk per the NCI criteria, and IKZF1del, all of which are markers of poor prognosis.

In conclusion, Pan-Cancer RNA sequencing identified all the different subgroups of BCP-ALL, except for high-hyperdiploid, hypodiploid, iAMP21, and whole gene (IKZF1, ERG, and PAX5) deletions. Previous associations between clinico-biological features and the different subgroups of ALL were mostly confirmed, and a gender bias toward males in PAX5-driven ALL and females in B-rest ALL was discovered.

The authors thank Maria Alessandra Salgado and Eduardo William Rodrigues Amâncio for clinical data management. The authors thank Lilian Girotto Zambaldi and Daniele Ribeiro Lucon for the excellent technical help with cytogenetics analysis, and Guilherme Giusti Navarro for helping with bioinformatics data management.

This work was supported by grants to J.A.Y. from Institute Ronald McDonald (project IRM 2016109) and CNPq (process 429811/2016-0). J.A.Y. received a Productivity fellowship from the Brazilian National Counsel of Technological and Scientific Development (CNPq, process 301596/2017-4 and 308399/2021-8). K.B.M. received financial support through FAPESP (grants 14/21704-9 and 14/50897-0). N.A.M., V.S.V., and G.L.C. received scholarships from Brazilian funding agencies (CAPES or CNPq).

Contribution: N.A.M. and J.A.Y. conceived and designed the study, performed the data interpretation, and wrote the manuscript; N.A.M. and P.Y.J. performed all targeted RNA-sequencing experiments; N.A.M. performed the validation assays by RT-PCR; P.Y.J., N.P.d.N., and G.L.C. contributed with GATA3 polymorphisms analysis by Sanger sequencing; N.A.M., V.S.V., and K.B.M. contributed with bioinformatics analysis; A.C.d.A. and S.R.B. were responsible for patients’ treatment and provided their demographic and clinical data; N.A.M. and J.A.Y. performed statistical analysis; and all authors revised the manuscript.

Conflict-of-interest disclosure: The authors declare no competing financial interests.

Correspondence: José Andrés Yunes, Centro Infantil Boldrini, Rua Gabriel Porto 1270, Campinas, SP 13083-210, Brazil; e-mail: andres@boldrini.org.br.