Transcriptomic classes of BCR-ABL1 lymphoblastic leukemia

Abstract

In BCR-ABL1 lymphoblastic leukemia, treatment heterogeneity to tyrosine kinase inhibitors (TKIs), especially in the absence of kinase domain mutations in BCR-ABL1, is poorly understood. Through deep molecular profiling, we uncovered three transcriptomic subtypes of BCR-ABL1 lymphoblastic leukemia, each representing a maturation arrest at a stage of B-cell progenitor differentiation. An earlier arrest was associated with lineage promiscuity, treatment refractoriness and poor patient outcomes. A later arrest was associated with lineage fidelity, durable leukemia remissions and improved patient outcomes. Each maturation arrest was marked by specific genomic events that control different transition points in B-cell development. Interestingly, these events were absent in BCR-ABL1⁺ preleukemic stem cells isolated from patients regardless of subtype, which supports that transcriptomic phenotypes are determined downstream of the leukemia-initialing event. Overall, our data indicate that treatment response and TKI efficacy are unexpected outcomes of the differentiation stage at which this leukemia transforms.

Fig. 1. Three molecular subtypes of BCR-ABL1 lymphoblastic leukemia.

a, Gene expression heatmap of three molecular subtypes in BCR-ABL1 lymphoblastic leukemia identified by consensus hierarchical clustering. Molecular subtypes, clustering confidence scores, BCR-ABL1 isoforms, clinical diagnoses and sample types are shown in tracks. A subset of hematopoietic lineage genes in each NMF component is shown. b, Heatmap of selected differentially expressed genes (FDR-adjusted P < 0.05 from pairwise comparison of subtypes) in stem/myeloid and B-lymphoid programs. Sample order and heatmap scale are the same as in a. Tracks on the left display from which pairwise comparison(s) each gene was derived. c, Proportions of blasts positive for myeloid (top) and B-lymphoid (bottom) lineage antigens by subtype. For cyMPO, the dashed line represents the diagnostic threshold of 10% and numbers in brackets represent outliers. Counts represent the number of primary leukemias assessed for the antigen. FDR-adjusted P values from Kruskal–Wallis tests are shown. In each boxplot in this paper, the middle line represents the median, the lower and upper edges of the box represent the first and third quartiles, the end of the lower whisker represents the smallest value at most 1.5× interquartile range from the lower edge of the box, the end of the upper whisker represents the largest value at most 1.5× interquartile range from the upper edge of the box and all data points were shown.

Fig. 2. Molecular subtypes are arrested at distinct stages of B-cell differentiation.

a, PCA of BCR-ABL1 lymphoblastic leukemias (circles) and normal hematopoietic cell types (triangles) using 2,000 most variable genes. b, Proportions of leukemias in each subtype with and without expression of clonally rearranged immunoglobulin light chain genes (IGK, IGL). P value is from Fisher’s exact test. c, Normalized counts of B-cell development genes, EBF1 and PAX5, by subtype (26 C1, 8 C2 and 23 C3 leukemias). d, Gene set score of B-cell receptor signaling by subtype (26 C1, 8 C2 and 23 C3 leukemias). e, Approximate positioning of the molecular subtypes in relation to B-cell differentiation stages. Transition from early pro-B to pro-B is accompanied by decreasing myeloid potential and increasing commitment to the B-lymphoid program. Timings of immunoglobulin heavy chain (IGH) and light chain (IGK, IGL) rearrangements are also shown. f, Single cells from scRNA-seq are annotated with their closest normal cell counterparts and visualized using UMAP. Exemplar cases from three subtypes are shown. Bars represent the proportions of cell type annotations. g, Proportions (median with interquartile range) of cell type annotations in each subtype (four Early-Pro, one Inter-Pro and four Late-Pro). Asterisks represent FDR-adjusted P < 0.05 from Wilcoxon rank-sum tests comparing Early-Pro and Late-Pro leukemias.

Fig. 3. Distinct cooperating genetic alterations define each molecular subtype.

a, Oncoprint of genetic alterations enriched in each subtype. Samples are in the same order as in Fig. 1a. Molecular subtype, sample type (diagnosis/relapse) and gene mutations are shown for each sample. Alteration frequencies in each subtype are shown on the right. IKZF1 deletions include monosomy 7. Copy number gains from trisomy 21 are not included in RUNX1 alteration frequencies. In hyperdiploid cases, trisomy 21 is a consequence of the hyperdiploid state. Fisher’s exact test was performed using only diagnostic samples or all samples. b, Recurrent IKZF1 deletions (blue lines) are generated by RAG-mediated recombination. RSS motifs near breakpoint clusters are shown. Right-side RSS are written in reverse complement and italicized. Bases that deviate from the canonical RSS (‘CACAGTG’) are in red. H3K4me3 ChIP–seq signal for GM12878 (B-lymphoblast cell line) from ENCODE is shown at the bottom. c, Proportions of BCR-ABL1-associated and transformation-related (Transform.) SVs with and without RSS motifs. d, Sequence logo of cryptic RSS heptamers from transformation-related SVs generated using WebLogo. e, Frequencies of leukemias with alterations in each gene. Colors correspond to 0, 1 or 2 RAG-mediated SVs. f, Distances between SV breakpoints and nearest RSS motifs (black) in transformation-related SVs form a negative binomial distribution (gray). g, Numbers of RAG-mediated recombinations in primary leukemias by deletion status of SLX4IP. Sig, statistical significance; Dx, diagnosis. ****P < 5 × 10⁻⁵ < ***P < 5 × 10⁻⁴ < **P < 0.005 < *P < 0.05.

Fig. 4. Origin of BCR-ABL1 lymphoblastic leukemia in a multipotent cell.

a, Simplified lineage hierarchy of B-lymphoblast and four colony types. Two potential scenarios for the origin of BCR-ABL1 are shown as follows: lymphoid progenitor origin (i) and stem cell origin (ii). Red outlines denote cell types that can harbor BCR-ABL1 under each scenario. b, Four patient samples with BCR-ABL1 positive colonies. For each colony type, numbers of BCR-ABL1 positive and negative colonies are shown. c, Timelines of clinical events in four patients with diagnosis and relapse leukemia pairs. Two circles correspond to the times of two leukemia biopsies, and circle colors represent molecular subtypes, Early-Pro (orange) and Late-Pro (blue). d, Phylogenetic relationships of four leukemia pairs based on point mutations. Lengths of trunk (shared) and branches (private) are proportional to number of mutations. Percentage of mutations in diagnostic leukemia that are shared by relapse leukemia is shown. The Jaccard index measures the degree of similarity between pairs. e, PCA of leukemic transcriptomes using NMF component genes. Diagnosis/relapse pairs are labeled and highlighted. f, PAX5 inactivation (box) and monosomy 7 are private to Ph53-D. Red and blue segments denote copy number gains and losses, respectively. g, RUNX1 mutation is private to Ph53-R. RUNX1 gene is on the minus strand and orange and blue represent G and C bases, respectively. A gray line represents a sequencing read, and a bar at the top represents the number of reads covering each position. Frequencies of reads with the c.445G>C mutation are shown. Lym, lymphoid progenitors; Mye, myeloid progenitors; GM, colony-forming unit (CFU)-granulocyte/macrophage; G, colony-forming unit-granulocyte; M, colony-forming unit-macrophage; E, burst-forming unit-erythroid; rel, relapse; D, diagnostic; R, relapse.

Fig. 5. Molecular subtypes predict treatment response and outcome.

a, Residual disease plot of log reduction in BCR-ABL1 transcript levels for each subtype. Each line represents a patient and each point represents a residual disease measurement. Dashed line denotes log reduction of 3, and log reduction of ≥3 is considered an MMR. Timings of first kinase domain mutation detections and BMTs are shown. b, First log reduction levels after induction therapy (at most 60 d from diagnosis) separated by subtype (20 Early-Pro, 12 Inter-Pro and 25 Late-Pro patients). Dashed line denotes log reduction of 3. P value from Kruskal–Wallis test is shown. c, Kaplan–Meier estimates of OS (top) and EFS (bottom) by subtypes. d, Kaplan–Meier estimates of OS by treatment type (imatinib-only versus switch to second- or third-generation TKI) for each subtype. e, Multivariable Cox proportional hazards model incorporating age (scaled by 10), sex, WBC (scaled by median), subtype, treatment type and interaction of subtype and treatment type (gray shade; n = 81). Colors represent variables associated with Inter-Pro (light blue) and Early-Pro (orange). Numbers on the left are P values for each variable. ND, not detected; ref, reference. **P < 0.005 < *P < 0.05.

Fig. 6. Kinase domain-independent mechanisms of treatment resistance.

a, Heatmap of average gene expression for cell quiescence and cell cycling gene sets. Sample order and subtype colors are the same as in Fig. 1a. b, Enrichment of cell quiescence gene set in the Early-Pro subtype (versus rest). c, Average log₂ expression (score) of G₁/S and G₂/M phase gene sets in single cells. Exemplar cases from three subtypes are shown. Dashed lines denote the G₁/S and G₂/M cutoffs (one standard deviation greater than mean). Bars represent proportions of cells in G₀, G₁/S and G₂/M phases. d, Proportion of cycling cells (non-G₀ phases) in each sample by subtype (four Early-Pro, one Inter-Pro and four Late-Pro leukemias). P value from the Wilcoxon rank-sum test comparing Early-Pro and Late-Pro samples is shown. e, Heatmap of average gene expression for various signaling gene sets. Sample order and subtype colors are the same as in Fig. 1a. f, Enrichment of BCR-ABL1 signaling and IL2/STAT5 signaling gene sets in the Early-Pro subtype (versus rest).

Fig. 7. Cell-of-transformation model for BCR-ABL1 lymphoblastic leukemia subtypes.

The disease originates from the rearrangement of BCR-ABL1 in a multipotent stem cell (HSC or MPP), which is the cell-of-origin. BCR-ABL1 increases cell proliferation and impairs DNA damage response. It then transforms by acquiring cooperating alterations in a downstream B-cell progenitor, which is the cell-of-transformation. Transformation-related events are mostly generated by RAG-mediated recombination. Inactivation of B-cell transcription factors and tumor suppressors causes a block in differentiation and evasion of programmed cell death, respectively. The molecular subtypes of BCR-ABL1 lymphoblastic leukemia arise from different target cells of transformation.

Extended Data Fig. 1. Flow sorting and sample clustering workflows.

a, Study design of the main cohort. b, Exemplar fluorescence-activated cell sorting scheme for leukemic blasts and T-cells. c, RNA-seq clustering approach incorporating non-negative matrix factorization (NMF) and consensus hierarchical clustering.

Extended Data Fig. 2. Three molecular subtypes of BCR-ABL1 lymphoblastic leukemia.

Consensus hierarchical clustering using 163 NMF component genes produced three molecular subtypes of BCR-ABL1 lymphoblastic leukemia. Rows represent the genes grouped into NMF components, and columns represent the 57 samples. Molecular subtypes, consensus clustering confidence scores, BCR-ABL1 isoforms, diagnoses (B-ALL, MPAL, CML-LBC), and sample types (primary/diagnostic, relapse) are shown in tracks.

Extended Data Fig. 3. Patient and sample characteristics of three transcriptomic subtypes.

Comparison of clinical and sample characteristics of three molecular subtypes. Plots are separated into main cohort (left; 23 C1, 8 C2, 22 C3), 3’-seq cohort (centre; 13 C1, 13 C2, 14 C3), and combined cohort (right; 36 C1, 21 C2, 36 C3). Counts represent numbers of patients/primary leukemias. Kruskal-Wallis test was used for comparing age at diagnosis and WBC, and Fisher’s exact test was used for the rest. PB, peripheral blood; BM, bone marrow; WBC, white blood cell.

Extended Data Fig. 4. Cell type annotation of single-cell RNA-seq data.

a, Single cells from scRNA-seq are annotated with their closest normal cell counterparts and visualized using UMAP. Six cases (3 Early-Pro and 3 Late-Pro) not shown in Fig. 2f are shown here. Bars represent proportions of annotated cell types. b, Proportions of annotated cell types in nine scRNA-seq samples. c, Shannon diversity index of cell type annotations in each single-cell RNA-seq sample. p-value from Wilcoxon rank-sum test comparing Early-Pro and Late-Pro subtypes is shown.

Extended Data Fig. 5. Frequent inactivation of lymphoid or B-cell development genes.

a, Oncoprint of alterations in lymphoid or B-cell development genes in decreasing order of recurrence rate. b, Histogram of the number of lymphoid or B-cell development genes altered in diagnostic leukemias. Mean number is 2.1. c, Frequency of homozygous losses of IKZF1 by subtype. p-value from Fisher’s exact test is shown. d, Frequency of dominant negative Ik6 deletion by subtype. p-value from Fisher’s exact test is shown. e, Gene expression score for stromal cell signature by subtype (26 Early-Pro, 8 Inter-Pro, 23 Late-Pro leukemias). f. Normalized counts (RNA-seq) of IL7R transcripts by subtype (26 Early-Pro, 8 Inter-Pro, 23 Late-Pro leukemias). g, Mutation loads (SNV and indel counts) by subtype. p-values from pairwise Wilcoxon rank-sum tests and Kruskal-Wallis test are shown. Data from primary leukemias with a median genome coverage of 15.8x are shown (12 Early-Pro, 3 Inter-Pro, 13 Late-Pro).

Extended Data Fig. 6. Detection and analysis of cryptic RSS sequences SV breakpoints.

a, Sequence logo of canonical RSS heptamers. b, Sequence logo of cryptic RSS heptamers from BCR-ABL1 SVs with RSS motifs (n = 20). Last 3 bases are not conserved for ‘GTG’. c, Distances between SV breakpoints and nearest RSS motifs for BCR-ABL1 SVs with RSS motifs. Unlike RSS motifs from transformation related SVs, they do not form a negative binomial distribution. d, Proportions of transformation SVs with and without RSS motifs by subtype.

Extended Data Fig. 7. Colony formation assay.

a, Schematic representation of the colony formation assay. b, Schematic representation of the nested PCR. Two breakpoints of the assayed SV are located >5 kb apart on the same chromosome or on different chromosomes (for example BCR in chr22 and ABL1 in chr9). PCR amplification is successful only when the assayed SV is present. c, Numbers of BCR-ABL1 positive and negative colonies in 14 samples. Subtype and BCR-ABL1 isoform (p210/p190) of each sample are shown.

Extended Data Fig. 8. Private and shared alterations in Ph53 at diagnosis and relapse.

a, Shared BCR-ABL1 translocation between Ph53-D and Ph53-R as evidenced by identical losses of flanking regions. b, SETD2 inactivation is private to Ph53-D. c, FIP1L1 deletion, gain of 18p, and trisomy 21 are private to Ph53-R. In A-C, blue represents copy number loss and red represents copy number gain relative to the reference T-cell genome. d, Nested PCR validation of SVs in Ph53. BCR-ABL1 translocation is shared, PAX5 and SETD2 inactivations are private to Ph53-D, and FIP1L1 deletion is private to Ph53-R. NTC, no template control. e, Sanger sequencing validation of C>G substitution at chr21:36252917 (RUNX1 c.445G>C, p.A149P) in Ph53-R (right). This mutation is not detected in Ph53-D (left).

Extended Data Fig. 9. Survival analysis by HBS1L status and by subtype-specific genomic alterations.

a, Survival analysis of patients grouped by HBS1L status (wildtype vs. deleted). b, Survival analysis of patients grouped by HBS1L status and molecular subtypes. Late-Pro and Inter-Pro patients are all HBS1L wildtype. c, Survival analysis of patients grouped by subtype-specific genomic alterations. Early-Pro alterations are HBS1L, RUNX1, EBF1, and/or Monosomy 7, and Inter-Pro alteration is two-copy loss of IKZF1. Kaplan-Meier estimates of overall survival (left) and event-free survival (right) are shown. Only the main cohort (n = 43 for survival analysis) is analyzed since these patients have corresponding genomic data. wt, wildtype; del, deleted.

Extended Data Fig. 10. Treatment resistance mechanisms in BCR-ABL1 lymphoblastic leukemia.

The degree of innate resistance to chemotherapy and TKI (that is imatinib) influences the probability of acquiring BCR-ABL1 KD mutations and relapsing. Late-Pro leukemias respond extremely well to induction and most patients achieve remission/MMR. Inter-Pro leukemias do not respond as well and most patients relapse without KD mutations. Early-Pro leukemias show minimal response to induction and most patients relapse with or without KD mutations. Cell quiescence, STAT5 signaling, and UPR signaling potentially contribute to treatment resistance in the Early-Pro and Inter-Pro subtypes.