Disease recurrence causes significant mortality in B-progenitor acute lymphoblastic leukemia (B-ALL). Genomic analysis of matched diagnosis and relapse samples shows relapse often arising from minor diagnosis subclones. However, why therapy eradicates some subclones while others survive and progress to relapse remains obscure. Elucidation of mechanisms underlying these differing fates requires functional analysis of isolated subclones. Here, large-scale limiting dilution xenografting of diagnosis and relapse samples, combined with targeted sequencing, identified and isolated minor diagnosis subclones that initiate an evolutionary trajectory toward relapse [termed diagnosis Relapse Initiating clones (dRI)]. Compared with other diagnosis subclones, dRIs were drug-tolerant with distinct engraftment and metabolic properties. Transcriptionally, dRIs displayed enrichment for chromatin remodeling, mitochondrial metabolism, proteostasis programs, and an increase in stemness pathways. The isolation and characterization of dRI subclones reveals new avenues for eradicating dRI cells by targeting their distinct metabolic and transcriptional pathways before further evolution renders them fully therapy-resistant.

Significance:

Isolation and characterization of subclones from diagnosis samples of patients with B-ALL who relapsed showed that relapse-fated subclones had increased drug tolerance and distinct metabolic and survival transcriptional programs compared with other diagnosis subclones. This study provides strategies to identify and target clinically relevant subclones before further evolution toward relapse.

See related video:https://vimeo.com/442838617

See related article by E. Waanders et al.

Despite significant advancements in the treatment of acute lymphoblastic leukemia (ALL), the disease recurs in 15% to 20% of pediatric and 40% to 75% of adult patients, with the prognosis for patients who relapse being dismal (1–3). Analysis of paired diagnosis and relapse ALL samples with increasingly broader and deeper genomic-sequencing methods has shown that classic Darwinian branching evolution of genomically distinct subclones is a hallmark of disease recurrence (4, 5). At both diagnosis and relapse, a single neoplasm may contain multiple genetic subclones related to each other through complex evolutionary trajectories (4–7). Although relapse may evolve from the predominant clone at diagnosis, in the majority of patients relapse arises from preexisting minor subclones within the diagnosis sample or from a rare ancestral clone (4–6). Although the population-level genetic analyses upon which these conclusions are based rely on computational inference of their evolutionary relationships from analysis of bulk leukemic cells, resolution of leukemic subclones at clonal levels, either through isolation of subclones in xenografts or from single-cell analysis, has largely substantiated these predictions (8, 9). Functional studies to explain therapy failure have mainly been undertaken by comparing diagnosis cells with those from relapse and identifying drug-resistance mechanisms present in relapse and absent at diagnosis. However, the properties of the relapse samples have been shaped by exposure to chemotherapy, causing further evolution and mutagenesis. Thus, two critical questions remain: What are the unique properties and mechanisms that contribute to the relapse fate of a particular diagnosis subclone prior to full evolution to drug-resistant relapse disease, and when does drug tolerance arise? Drug tolerance may arise stochastically through genetic or epigenetic mechanisms prior to exposure to therapy, and be selected for by both cell-autonomous and non–cell autonomous processes (10–13). Alternatively, therapy may induce genomic aberrations that are then selected for during disease progression, particularly if such alterations reduce leukemic cell fitness (14) during disease establishment and prior to the administration of therapy (14, 15, 16). Without isolation of the subclones that contribute to disease progression from diagnosis samples, it is not possible to answer these questions and uncover the cellular and molecular properties that explain their differing subclonal fates and drug tolerance.

Many therapy resistance mechanisms have been implicated in B-progenitor ALL (B-ALL), including acquisition of stemness programs, dormancy, the protective role of the niche, and the acquisition of resistance driver mutations (14, 17–21). However, the ability to evade drug treatment is only one prerequisite of relapse; surviving cells must also possess significant clonal regenerative capacity to regrow or reproduce disease. For many human cancers, only rare fractions of malignant cells are capable of significant clonal propagation as detected by xenograft-based cancer or leukemia-initiating cell (L-IC) assays (10). Indeed, methods to propagate primary human leukemia samples were first undertaken in B-ALL with patient-derived xenografts (PDX) recapitulating many features of the patient's disease (22, 23). Subsequent studies used a limiting dilution approach to generate xenografts, thereby tracking the growth properties and genetic identity of single L-ICs (9). Thus, depending on the cell dose transplanted, rare subclones with poorer competitive repopulation properties, compared with more aggressive or numerous L-ICs within the sample, could be identified. By genetic analysis of such xenografts, evidence was found for the existence of ancestral and/or minor subclones, proving that branching evolution and clonal diversity occur at the level of L-ICs; however, these studies did not test for relapse-fated subclones (8, 9). Conceptually, L-ICs with the capacity for clonal propagation should serve as the units of selection during disease progression because only mutations accumulating in clonal propagating cells are relevant for further disease evolution. Thus, the relationship between disease progression and properties of stemness is an active area of investigation (10, 24). The pairing of xenografting assays with genomic studies provides a unique opportunity to enrich for the cellular reservoirs of relapse. However, a direct test of this concept through paired diagnosis and relapse xenografting in B-ALL has been limited (23, 25, 26). Paired studies in acute myeloid leukemia (AML) and T-lineage acute lymphoblastic leukemia (T-ALL) have demonstrated that xenografts can capture distinct subclones present within the diagnosis samples of relapsing patients, some genetically more closely related to the predominant diagnosis subclone and others to a corresponding relapse sample (27–30). These results suggest that latent relapse-initiating subclones possess competitive growth advantages when assessed by xenografting, but remain suppressed by the predominant clone in the patient, making them difficult to study.

Here, we undertook a combined functional and genomic analysis of 14 genetically distinct paired diagnosis and relapse B-ALL patient samples to isolate latent subclones within the diagnosis sample that initiate relapse. Isolation of subclones enabled functional analysis of their growth and drug response properties as well as molecular analysis of their transcriptomic profiles.

Isolation of Relapse-Initiating Subclones in B-ALL

Whole-exome sequencing (WES) and SNP microarray analysis of 6 adult patients and 8 pediatric patients with B-ALL of varying genetic subtypes were undertaken to identify somatic single-nucleotide variants (SNV), insertion–deletion mutations (indels), and DNA copy-number alterations (CNA). The patients encompassed a range of cytogenetic subtypes and varied in the length of their disease remission (range 5.88–94.8 months; Supplementary Table S1). Diagnosis samples had a median of 24 somatic SNV/indels (range 7–100) and 13.5 CNA (range 1–51), whereas relapse samples contained a median of 39.5 SNV/indels (range 22–405) and 16.5 CNA (range 2–58; Supplementary Table S1). Leukemic variants were confirmed by targeted sequencing using a custom capture array of the variants identified by WES (Fig. 1A). Targeted sequencing and WES data were merged, resulting in a coverage of approximately 350×. Computational analysis of the variant allele frequencies (VAF) of leukemic variants comparing diagnosis and relapse samples predicted that the origin of relapse arose from a minor subclone present at diagnosis in 10 patients (patients 1, 3–7, 9, and 12–14; Fig. 1B; Supplementary Table S1) and through further evolution of the major diagnostic subclone in 4 patients (patients 2, 8, 10, and 11; Supplementary Table S1). Broadly, these findings of evolutionary origins are representative of the much larger analysis of 92 paired samples that reports the mutational landscape and patterns of clonal evolution of relapsed childhood ALL as described in Waanders and colleagues (7).

Figure 1.

PDXs capture clonal diversity present in paired diagnosis and relapse B-ALL samples. A, Experimental schematic. PDXs were generated for 6 adult and 8 pediatric B-ALL patient samples at diagnosis and relapse by intrafemoral transplantation of sorted leukemic blasts into 30 irradiated NSG mice in a limiting dilution assay. Mice were sacrificed 20 to 30 weeks post-transplant and their engraftment was assessed by flow cytometry. Patient samples were also subjected to genomic analysis including whole-exome sequencing (WES). Variants identified from the WES of patient samples at either time point were used to create custom capture baits for targeted sequencing at a deeper depth in the patient samples and their corresponding PDX. PDXs representing varying clones were identified. B, Schematic representation of the results obtained by mutational clustering of variants based on the variant allele frequencies (VAF) at diagnosis (x-axis) and relapse (y-axis) of patient 1 in 2-D VAF plots showing evolution from a minor subclone as depicted. Each dot represents a variant. Shared variants are shown in gray clusters (clusters a and c). Diagnosis- and relapse-specific variants are shown in the blue cluster (cluster b) and red clusters (clusters d and e), respectively. C, Heat map of the VAFs of leukemic variants at diagnosis, relapse, and in their corresponding PDXs for patient 1. Variant classes are labeled with their class and a letter corresponding to the clusters illustrated in B.

Figure 1.

PDXs capture clonal diversity present in paired diagnosis and relapse B-ALL samples. A, Experimental schematic. PDXs were generated for 6 adult and 8 pediatric B-ALL patient samples at diagnosis and relapse by intrafemoral transplantation of sorted leukemic blasts into 30 irradiated NSG mice in a limiting dilution assay. Mice were sacrificed 20 to 30 weeks post-transplant and their engraftment was assessed by flow cytometry. Patient samples were also subjected to genomic analysis including whole-exome sequencing (WES). Variants identified from the WES of patient samples at either time point were used to create custom capture baits for targeted sequencing at a deeper depth in the patient samples and their corresponding PDX. PDXs representing varying clones were identified. B, Schematic representation of the results obtained by mutational clustering of variants based on the variant allele frequencies (VAF) at diagnosis (x-axis) and relapse (y-axis) of patient 1 in 2-D VAF plots showing evolution from a minor subclone as depicted. Each dot represents a variant. Shared variants are shown in gray clusters (clusters a and c). Diagnosis- and relapse-specific variants are shown in the blue cluster (cluster b) and red clusters (clusters d and e), respectively. C, Heat map of the VAFs of leukemic variants at diagnosis, relapse, and in their corresponding PDXs for patient 1. Variant classes are labeled with their class and a letter corresponding to the clusters illustrated in B.

Close modal

To gain insight into the genetic diversity at the level of L-ICs and uncover rare and/or outcompeted clones, purified leukemic blasts from primary diagnosis and relapse samples were transplanted intrafemorally in a limiting dilution assay into 872 NOD.CB17-PrkdcscidIl2rgtm1Wjl/Szj (NSG) mice to generate primary PDXs. The frequency of L-ICs varied widely between samples, with L-IC ranges corresponding to those described previously (31, 32); however, paired analysis between diagnosis and relapse samples did not show a consistent trend in L-IC enrichment at either time point (Supplementary Table S2). This collection of PDXs that were engrafted [n = 402, average 28 per patient (range 11–58), 13 per sample (range 0–27)] was used to further study the clonal landscape present in the patient samples. To determine whether PDXs captured the clonal diversity and disease evolution present in the patient samples, human cells were isolated from the bone marrow (BM), spleen, and central nervous system (CNS) of engrafted mice and subjected to targeted sequencing using the custom capture array designed for the patient samples (Fig. 1A). PDXs with sufficient human engraftment for genotyping (>10% human chimerism, 372 PDXs total of 402 engrafted PDXs, average 26 PDXs per patient) were analyzed. Individual PDXs from the same patient sample were often found to vary in their clonal composition in terms of the presence and frequency of variants, suggesting that the L-ICs initiating the grafts derived from genetically diverse subclones often capture the totality of the clonal diversity present in the diagnosis patient sample (Fig. 1C).

Leukemic variants were classified on the basis of the VAF of the variants in the bulk patient diagnosis and relapse samples from which the PDXs were generated: preserved variants (VAF > 30% in both diagnosis and relapse samples, or preserved between samples); diagnosis-specific variants (present at diagnosis and absent at relapse); latent variants (present at diagnosis with VAF < 30% and increasing at relapse); relapse-specific variants (absent at diagnosis and present at relapse; Supplementary Table S3). The limit of detection of our combined sequencing was a VAF of 1%. This analysis revealed three patterns of engraftment in PDXs derived from diagnosis samples (termed dPDXs). In 10 of 13 engrafting diagnosis samples (76.9%, patients 1–7, 9, 12, and 14), latent variants were enriched in dPDXs (>10% increase in VAF in dPDXs) as compared with the diagnosis patient sample demonstrating the regenerative potential of clones marked by these variants (Fig. 2A; Supplementary Fig. S1A and S1B). In 4 of 13 patients (patients 5, 6, 11, and 13), dPDXs were generated where leukemia cells bore relapse-specific variants establishing their existence within the diagnosis sample, despite being at levels below the detectable limit of the sequencing in the patient sample (Fig. 2B). Finally, in one patient (patient 8), whose relapse was predicted to evolve from the major diagnosis clone based on analysis of the bulk patient samples, only diagnosis clones engrafted in dPDXs (Supplementary Fig. S1C). This patient carried an ETV6–RUNX1 translocation and had the longest remission, relapsing 8 years after the initial B-ALL diagnosis. Our approach of generating PDXs with differing cell doses was instrumental in showing the existence of subclones at diagnosis that bear latent (patient 2) or relapse-specific variants indicative of ancestral clones (patient 11) in 2 patients in which genetic analysis of bulk diagnosis and relapse samples had predicted evolution from the major diagnosis clone (Fig. 2B; Supplementary Table S3). Only one diagnosis patient sample (patient 10), a sample predicted to arise from evolution of the major diagnostic clone, did not engraft at the highest transplanted cell dose (250,000 cells). Therefore, xenografting added considerable new insight into the subclonal repertoire of L-ICs in these patients and their evolutionary fates and patterns.

Figure 2.

PDXs enrich for latent diagnosis clones. A and B, Heat maps of VAF of the SNV and indel leukemic variants identified by whole-exome and targeted sequencing in diagnosis/relapse patient samples and PDXs, respectively. Variants are clustered as preserved (present in diagnosis and relapse patient samples), diagnosis specific (present in diagnosis patient sample and absent in relapse patient sample), latent (present in diagnosis patient sample with VAF < 0.3 and expanding in relapse sample), and relapse specific (present in relapse patient sample and absent in diagnosis patient sample). PDXs are ordered in decreasing number of transplanted cell doses. A, Representative heat map for the selection of latent variants in diagnosis PDX as observed in patient 9. B, Representative heat map of patient 11 displaying the identification of a relapse-specific variant undetectable in the patient diagnosis sample but present in diagnosis PDX.

Figure 2.

PDXs enrich for latent diagnosis clones. A and B, Heat maps of VAF of the SNV and indel leukemic variants identified by whole-exome and targeted sequencing in diagnosis/relapse patient samples and PDXs, respectively. Variants are clustered as preserved (present in diagnosis and relapse patient samples), diagnosis specific (present in diagnosis patient sample and absent in relapse patient sample), latent (present in diagnosis patient sample with VAF < 0.3 and expanding in relapse sample), and relapse specific (present in relapse patient sample and absent in diagnosis patient sample). PDXs are ordered in decreasing number of transplanted cell doses. A, Representative heat map for the selection of latent variants in diagnosis PDX as observed in patient 9. B, Representative heat map of patient 11 displaying the identification of a relapse-specific variant undetectable in the patient diagnosis sample but present in diagnosis PDX.

Close modal

Genetic Analysis of Xenografts Provides Insight into the Dynamics of Subclone Evolution

To gain insight into the evolutionary relationships and processes underlying the divergence of subclones, we undertook two approaches: population phylogenetic analysis to examine the genetic similarity between xenografts and patient samples; and generation of mutational trees to reconstruct the clonal hierarchies. We began by using a population genetics approach where phylogenetic analysis was performed for the patient samples and the clones isolated from their PDXs to depict the evolutionary relationships between xenografts and the patient samples (Fig. 3A and B; Supplementary Fig. S1D–S1G). These analyses demonstrate the ability of PDXs to capture diagnosis subclones more closely related to the diagnosis patient sample, as well as subclones present in the diagnosis patient sample that are on the trajectory to relapse. We term these partially evolved subclones that go on to cause relapse diagnosis Relapse Initiating clones (dRI). dRIs were identified in all engrafting diagnosis patient samples except for the patient with the longest time of remission (patient 8). Thus, our xenografting approach enabled the capture and isolation of subclones, including those on the evolutionary trajectory of relapse, allowing for further functional analysis and aiding in the determination of mutational acquisition.

Figure 3.

PDXs identify relapse-fated clones in diagnosis patient samples. Phylogenetic analysis showing relationship of patient 9 (A) and patient 11 (B) patient samples and xenografts, based on VAF of leukemic variants. The distances between symbols on the tree were estimated by a nearest neighbor joining method and represent the degrees of relation between them (Minkowski's distance). Circles represent patient samples and triangles represent PDX; blue represents diagnosis and red represents relapse. Diagnosis clones on the trajectory to relapse were termed dRI and are indicated by a box with a hatched border.

Figure 3.

PDXs identify relapse-fated clones in diagnosis patient samples. Phylogenetic analysis showing relationship of patient 9 (A) and patient 11 (B) patient samples and xenografts, based on VAF of leukemic variants. The distances between symbols on the tree were estimated by a nearest neighbor joining method and represent the degrees of relation between them (Minkowski's distance). Circles represent patient samples and triangles represent PDX; blue represents diagnosis and red represents relapse. Diagnosis clones on the trajectory to relapse were termed dRI and are indicated by a box with a hatched border.

Close modal

For our second approach, we used the extensive mutational information garnered from the 372 PDXs to supplement the patient analysis. For each patient, computational analysis using PairTree (Wintersinger and colleagues, in preparation) was undertaken to group leukemic variants with similar VAFs across numerous PDX and bulk samples to define genetically unique mutational populations. We then constructed mutational trees to describe the evolutionary history of each leukemia and the order of mutational acquisition (Fig. 4; Supplementary Figs. S2A–S2C, S3A–S3C, S4A–S4E, and S5A–S5C; Supplementary Table S3). In 7 of 14 patients, the PDXs helped separate the mutational acquisition into multiple distinct populations that could not be resolved using the sequencing data of the primary patient samples alone. For instance, in patient 1, analysis of the xenografts revealed the existence of four distinct subclonal populations (populations 10, 12, 14, and 15) at relapse that, when only the patient samples were considered, had been predicted to be only a single subclone (population G; Fig. 4AD). This analysis further enabled comparisons of the state of the evolutionary trajectory of the subclones captured within each individual PDX (denoted by a numbered PDX referring to a specific mouse) to the patient samples from which they were derived, such as dPDX7 in patient 1 whose clonal composition is more closely related to the relapse than the diagnosis sample (Fig. 4C and D). In addition, the PDXs aided in subdividing single mutational populations into more than a single linearly related population or branch in 11 of 14 patients. For example, population G in patient 9 was expanded into two separate populations (populations 8 and 9) when PDXs were included (Fig. 4E and F; Supplementary Figs. S2A–S2C, S3A–S3C, S4A–S4E, and S5A–S5C).

Figure 4.

Generation of mutational trees from the combined genomic data of xenografts and bulk diagnosis and relapse patient samples. Mutational trees of variants clustered to form populations using the PairTree algorithm. Nodes in mutational trees are divided in half, with the intensity of blue in the left half indicating the frequency of the population's variants at diagnosis, and the intensity of red in the right half showing the frequency of the population's variants at relapse. Color intensity shows subclonal prevalence as noted in the legend of A and applicable to all other trees except C. A, Mutational tree of patient 1 derived by analysis of the patient samples (diagnosis and relapse) alone. Mutational populations identified from bulk patient samples alone are denoted by a square node labeled with an alphabetical letter. B, Combined mutational tree derived from the variant analysis of both patient 1 patient samples and all their generated xenografts. Mutational populations derived from combined patient and xenograft analysis are represented by circular nodes labeled with numerals. The mutational populations identified using the patient samples alone in A are overlaid on the tree as boxes labeled with their corresponding alphabetical letter. This identifies instances where single populations in A correspond to multiple populations when xenografts are included (i.e., mutational population G). C, Combined mutational tree of patient 1 shaded to indicate the prevalence of variants in dPDX 7 (instead of diagnosis and relapse) demonstrating that this PDX is composed primarily of variants of the relapse lineage. D, Presence of identified mutational populations in patient samples and representative xenografts. Mutational populations (Pop.) are displayed on the y-axis and individual patient samples or xenografts are displayed on the x-axis. The height of the population bar represents the prevalence of the lineage in the patient sample (Pt.) or xenograft. E, Mutational tree, similar to A, derived from patient samples alone of patient 9. F, Combined mutational tree, similar to B, derived jointly from patient 9 patient samples and xenografts. Subclonal prevalences of populations 2 to 5 are shown, indicating the absence of diagnosis populations 2 to 4 and presence of population 5 (the first node in the relapse lineage branch) in all dPDX. G, Mutational tree, similar to A, of patient 11 derived from patient samples alone. H, Combined mutational tree, similar to B, of patient 11 derived from patient samples and xenografts. The prevalence of mutational population 3 is displayed, highlighting its absence in the diagnosis patient and its detection in only two dPDXs.

Figure 4.

Generation of mutational trees from the combined genomic data of xenografts and bulk diagnosis and relapse patient samples. Mutational trees of variants clustered to form populations using the PairTree algorithm. Nodes in mutational trees are divided in half, with the intensity of blue in the left half indicating the frequency of the population's variants at diagnosis, and the intensity of red in the right half showing the frequency of the population's variants at relapse. Color intensity shows subclonal prevalence as noted in the legend of A and applicable to all other trees except C. A, Mutational tree of patient 1 derived by analysis of the patient samples (diagnosis and relapse) alone. Mutational populations identified from bulk patient samples alone are denoted by a square node labeled with an alphabetical letter. B, Combined mutational tree derived from the variant analysis of both patient 1 patient samples and all their generated xenografts. Mutational populations derived from combined patient and xenograft analysis are represented by circular nodes labeled with numerals. The mutational populations identified using the patient samples alone in A are overlaid on the tree as boxes labeled with their corresponding alphabetical letter. This identifies instances where single populations in A correspond to multiple populations when xenografts are included (i.e., mutational population G). C, Combined mutational tree of patient 1 shaded to indicate the prevalence of variants in dPDX 7 (instead of diagnosis and relapse) demonstrating that this PDX is composed primarily of variants of the relapse lineage. D, Presence of identified mutational populations in patient samples and representative xenografts. Mutational populations (Pop.) are displayed on the y-axis and individual patient samples or xenografts are displayed on the x-axis. The height of the population bar represents the prevalence of the lineage in the patient sample (Pt.) or xenograft. E, Mutational tree, similar to A, derived from patient samples alone of patient 9. F, Combined mutational tree, similar to B, derived jointly from patient 9 patient samples and xenografts. Subclonal prevalences of populations 2 to 5 are shown, indicating the absence of diagnosis populations 2 to 4 and presence of population 5 (the first node in the relapse lineage branch) in all dPDX. G, Mutational tree, similar to A, of patient 11 derived from patient samples alone. H, Combined mutational tree, similar to B, of patient 11 derived from patient samples and xenografts. The prevalence of mutational population 3 is displayed, highlighting its absence in the diagnosis patient and its detection in only two dPDXs.

Close modal

Genetically Diverse Subclones Have Differing Xenograft Repopulation Kinetics

The mutational analysis of the xenografts performed using PairTree helped illuminate the competitive differences in xenograft repopulation kinetics of specific subclones. We compared the predominant diagnosis mutational populations in the patient sample to the mutational populations that engrafted in the xenografts. In all but four patient samples, the dPDX captured the diversity of clones present in the diagnosis samples (Fig. 4; Supplementary Figs. S2A–S2C, S3A–S3C, S4A–S4E, and S5A–S5C). In 2 of these 4 patients (patients 7 and 11), minor diagnosis-specific mutational populations, including a population harboring missense NRAS mutations in patient 7, did not engraft in the dPDX mice (Supplementary Fig. S4A; Supplementary Table S3). In the other two patients (patients 9 and 12), all the dPDXs were initiated from a minor population in the diagnosis sample ancestral to the relapse (Fig. 4E and F; Supplementary Figs. S4C and S5A). These data show the enhanced competitive repopulation properties of these dRI clones. In patient 9, all dPDXs were initiated from mutational population 5 that corresponded to only 22% of the diagnosis leukemic cells (initially described as population E in the patient sample alone mutational tree) in comparison with the mutational population 2 lineage that represented 71% of the diagnosis leukemic cells (initially described in population B; Fig. 4E and F). In the dPDXs, population 5 displayed a selective advantage representing at least 93% of the leukemic cells, outcompeting the population 2 lineage (encompassing populations 2–4) and recapitulating the evolutionary dynamics of the patient relapse, where 98% of the leukemic cells came from the population 5 lineage (Fig. 4F). In contrast, in 5 patients (patients 1, 4, 6, 7, and 11), the relapse-fated mutational population represented only a small percentage of the engrafting cells in most dPDXs and only rose to predominance in rare mice. For example, in patient 11, relapse-fated mutational population 3, initially described in population C, corresponded to only 1% of leukemic cells at diagnosis and remained at similarly low levels in most dPDXs. However, in two dPDXs transplanted with a low, near limiting cell dose (10,000 cells), population 3 rose to 92% and 74% of leukemic cells, distinguishing it from population 4 and highlighting the importance of the limiting dilution approach to capture rare subclones by isolating them from other subclones with higher competitive capacity (Fig. 4G and H; Supplementary Fig. S4E). Population 3 was defined by a single variant in the 3′UTR of AHNAK, a gene whose expression has previously been implicated in disease relapse in T-ALL (33). Thus, the xenografting strategy permitted integration of functional and genomic information, providing insight into the distinct competitive growth properties of functionally defined L-IC subclones. Taken together, the phylogenetic analysis and mutational evolutionary trees provide strong evidence that relapse-fated subclones were already present in the diagnosis sample, confirming prior studies of several human acute leukemic diseases and in silico predictions (4, 5, 28, 29).

Genetically Diverse Subclones Have Differing Immunophenotypic and Migratory Properties

To further examine whether genetically distinct subclones also possessed variation in their immunophenotype or functional properties that might explain why some are fated to contribute to relapse and others are not, we examined their differentiation, growth, and migration properties. We first interrogated the differential properties of the diagnosis clones isolated from patient 9 in which the presence of latent variants, including the known relapse driver CREBBP (21), was segregated between clones. Differences in immunophenotype were observed for dRI-PDXs enriched for these latent variants that were present at low VAFs in the diagnosis patient sample (1%–11%; Fig. 5A and B). The dRI-PDXs exhibited a CD45CD34+ immunophenotype reflective of the relapse sample, whereas others contained both a CD45CD34+ and a CD45dimCD34+ population as observed in the diagnosis patient sample (Fig. 5A). Relapse PDX (rPDX), like the relapse sample itself, were CD45 (Fig. 5A). dPDXs that were engrafted primarily with CD45dim or both immunophenotypic populations had a genotype that was more closely related to the diagnosis sample (Fig. 5B). The differences noted in immunophenotype appeared to reflect the rise of the relapse-fated clone and aided in identifying different subclones present within a single xenograft. To validate that the difference in immunophenotype segregated the predominant diagnosis clone from the minor subclone that seeded the relapse, we isolated the two immunophenotypic populations by flow-sorting cells from 6 dPDXs and subjected the populations to targeted sequencing using the custom capture array designed for the patient samples. This analysis confirmed the enrichment of the dRI clone in the CD45 population, whereas latent variants were absent or very rare (<6%) in the CD45dimCD34+ population (Fig. 5B). Of note, change in CD45 expression between diagnosis and relapse timepoints was identified in 5 of our patient samples. PDX from patient 9 showed markedly reduced leukemic dissemination from the injected femur to other hematopoietic sites (Supplementary Fig. S6A).

Figure 5.

Competitive dRI-PDX clones identified in diagnosis PDXs. A, Flow cytometry analysis of patient samples and representative dPDX and rPDX of patient 9 display the presence of two different immunophenotypic populations: a CD45dimCD34+ and a CD45CD34+. B, Targeted sequencing of the dPDX revealed variability in the VAF of latent variants that corresponded with the shift in immunophenotypic populations. Cell sorting for immunophenotypic populations followed by targeted sequencing revealed the isolation of the latent variants in the CD45CD34+ population. C, RNA-seq analysis of differentially expressed genes (adjusted p-value of <0.05 and absolute log2 fold change of >1) between the two dPDX (dRI-PDX clone 2; dRI-PDX clone 1CREBBP_WT) versus rPDX for Patient 9. Relative expression was generated from variance stabilized normalized counts. D, Enrichment map of gene sets differentially enriched in patient 9 dRI-PDX clone 1CREBBP_WT versus dRI-PDX clone 2 (FDR q value ≤ 0.05) and dRI-PDX clone 1CREBBP_WT versus rPDX clone 3 (FDR q value ≤ 0.05). Node size is proportional to the number of genes included in each gene set (minimum 10 genes/gene set). Gray and red edges indicate gene overlap. Green node: enrichment in dRI-PDX clone 1CREBBP_WT (positive NES). Purple node: enrichment in dRI-PDX clone 2 and/or rPDX clone 3 (negative NES). Autoannotate app in Cytoscape was used to automatically annotate clusters (black squares). NES, normalized enrichment score; DDR, DNA damage response; TCR, T-cell receptor; UPR, unfolded protein response; cGMP, cyclic guanosine monophosphate; DCC, deleted in colorectal cancer gene; NO, nitric oxide; CSK, C-terminal Src kinase; CFTR, cystic fibrosis transmembrane conductance regulator; SPRY, Sprouty gene. E, Human purified cells from primary dPDX and dRI-PDX from patient 4 were transplanted into secondary NSG recipients. Mice were monitored for peripheral blood human chimerism until mean blood levels reached greater than 10%, revealing different kinetics of chimerism between dPDX 7 and dRI-PDX 11. Symbols represent the mean chimerism value of dPDX 7 (n = 20 mice) and dRI-PDX 11 (n = 16 mice) and bars represent SD.

Figure 5.

Competitive dRI-PDX clones identified in diagnosis PDXs. A, Flow cytometry analysis of patient samples and representative dPDX and rPDX of patient 9 display the presence of two different immunophenotypic populations: a CD45dimCD34+ and a CD45CD34+. B, Targeted sequencing of the dPDX revealed variability in the VAF of latent variants that corresponded with the shift in immunophenotypic populations. Cell sorting for immunophenotypic populations followed by targeted sequencing revealed the isolation of the latent variants in the CD45CD34+ population. C, RNA-seq analysis of differentially expressed genes (adjusted p-value of <0.05 and absolute log2 fold change of >1) between the two dPDX (dRI-PDX clone 2; dRI-PDX clone 1CREBBP_WT) versus rPDX for Patient 9. Relative expression was generated from variance stabilized normalized counts. D, Enrichment map of gene sets differentially enriched in patient 9 dRI-PDX clone 1CREBBP_WT versus dRI-PDX clone 2 (FDR q value ≤ 0.05) and dRI-PDX clone 1CREBBP_WT versus rPDX clone 3 (FDR q value ≤ 0.05). Node size is proportional to the number of genes included in each gene set (minimum 10 genes/gene set). Gray and red edges indicate gene overlap. Green node: enrichment in dRI-PDX clone 1CREBBP_WT (positive NES). Purple node: enrichment in dRI-PDX clone 2 and/or rPDX clone 3 (negative NES). Autoannotate app in Cytoscape was used to automatically annotate clusters (black squares). NES, normalized enrichment score; DDR, DNA damage response; TCR, T-cell receptor; UPR, unfolded protein response; cGMP, cyclic guanosine monophosphate; DCC, deleted in colorectal cancer gene; NO, nitric oxide; CSK, C-terminal Src kinase; CFTR, cystic fibrosis transmembrane conductance regulator; SPRY, Sprouty gene. E, Human purified cells from primary dPDX and dRI-PDX from patient 4 were transplanted into secondary NSG recipients. Mice were monitored for peripheral blood human chimerism until mean blood levels reached greater than 10%, revealing different kinetics of chimerism between dPDX 7 and dRI-PDX 11. Symbols represent the mean chimerism value of dPDX 7 (n = 20 mice) and dRI-PDX 11 (n = 16 mice) and bars represent SD.

Close modal

The ability to use differences in immunophenotype to segregate clones enabled us to select evolutionarily related dRI subclones for RNA-sequencing (RNA-seq) analysis from patient 9. We compared the expression profiles of two dRI clones, an ancestral and a daughter clone (clones 1 and 2, respectively), to a representative relapse clone (clone 3). The dRI clone most genetically similar to the relapse, bearing latent leukemic variants including 3 variants in CREBBP (2 missense and 1 silent) as well as variants in TCS2, NRAP, PLXNA4, and PRINS (dRI-PDX clone 2) had an expression profile very closely related to representative relapse clones with only 24 differentially expressed genes (Fig. 5C and D; Supplementary Table S4). In comparison, the ancestral clone not harboring CREBBP (dRI-PDX clone 1CREBBP_WT) showed 479 differentially expressed genes compared with the relapse clones, suggesting that the majority of transcriptional changes that are seen at relapse had already occurred in dRI-PDX clone 2 at diagnosis prior to exposure to chemotherapy (Supplementary Table S4). Broadly, the changes in gene expression of dRI-PDX clone 2 in relation to dRI-PDX clone 1CREBBP_WT centered on dysregulation of histone variants and inflammation-related genes as well as downregulation of genes involved in morphogenesis (TGFβ signaling, GATA3), the T-cell receptor and B-cell receptor pathways (CD19, CD79A, CD4, etc.) and antigen processing and presentation pathway, while heatshock response (HSR) and unfolded protein respone (UPR) pathways were enriched (Supplementary Table S4). Thus, the early acquisition of additional leukemic variants in this relapse-fated subclone, including the CREBBP mutations, caused significant changes in gene expression. Interestingly, tyrosine kinases such as EPHB2, LTK, and ERBB2 were upregulated after acquisition of CREBBP variants in both dRI-PDX clone 2 and the relapse clone, suggesting possible vulnerabilities to tyrosine kinase inhibitors. Despite this striking change, the dRI subclone remained dormant within the patient for many years, as this patient displayed a long remission (4 years).

PDXs from the majority of samples showed extensive migration of leukemic cells to other hematopoietic sites and other tissues including the spleen, CNS, and peripheral blood (Supplementary Fig. S6B and S6C). Because patients with B-ALL may present with leukemic dissemination to the CNS and testes where they can provide a sanctuary for relapse, we investigated whether there were differences in the dissemination of subclonal populations to each site. Targeted sequencing analysis showed that there was genetic discordance between CNS and bone marrow in 40% (44 of 111) of xenografts and between spleen and BM in 17.8% (28 of 157) of xenografts (Supplementary Table S5). In one patient (patient 7), our analysis of the genetic discordance revealed the presence of a dRI clone engrafting in the CNS (VAF > 39%) of a xenograft transplanted with a high cell dose (250,000 cells), in which the clone was a minor clone barely detectable in the BM (VAF of < 3%) and outcompeted by a more predominant diagnosis clone (Supplementary Fig. S6B). The identification of the dRI subclone in the CNS of dPDX is consistent with the ability of relapse-fated cells to disseminate and cause disease recurrence in the CNS of patients with B-ALL. We also noted the occurrence of a difference in peripheral blood dissemination of dRI in an additional patient (patient 4), with delayed mobilization of the dRI-PDX subclone as compared with the representative diagnosis clone (Fig. 5E). Collectively, the combined functional and genotyping studies show that individual subclones possess distinct immunophenotypic, competitive repopulation, proliferative, and migration properties and that dRI subclones can already possess distinct biological properties at diagnosis, before exposure to chemotherapeutic agents.

dRI Subclones Display Differential Response to Chemotherapeutic Agents

To directly test the functional properties of dRI subclones for their ability to survive and contribute to relapse, we compared the drug sensitivities of individual subclones for 5 of the genetically distinct patients. Multiple secondary PDXs were generated from dPDXs with known dRI clones (dRI-PDX), predominant diagnosis clones (dPDX), or representative relapse clones (rPDX). Treatment of patients with B-ALL includes combination chemotherapy with supportive care; however, it is not possible to replicate human therapy precisely in xenografts. Our interest was to determine whether there was any variation in the responsiveness of different subclones to individual drugs used in these chemotherapeutic protocols. Following engraftment, PDXs received single-agent treatments of dexamethasone, vincristine, l-asparaginase, or saline for 4 weeks (Fig. 6A and B; Supplementary Fig. S7A–S7F). Differences in therapeutic responses to one or more drugs between dRI and representative diagnosis clones were observed for four patient samples. In three patient samples (patient 1, 6, 7), dPDXs harboring a dRI clone (dRI-PDX) were compared with dPDXs repopulated with the predominant diagnosis clone and demonstrated reduced sensitivity to 2 or 3 of the 3 chemotherapy agents tested (Fig. 6A and B; Supplementary Fig. S7A–S7D). Reduced sensitivity to a single chemotherapeutic agent was also observed in one additional patient (patient 4; dexamethasone, significant in injected femur and trend in BM and spleen; Supplementary Fig. S7E). In contrast, there was no difference in the therapeutic sensitivity of two dRI-PDXs from patient 11 defined by the presence of the AHNAK 3′UTR variant. This patient has a strong driver KMT2A (MLL) translocation, and the presence of the single relapse-specific variant does not appear to confer sufficient evolution of the leukemia to alter the therapeutic sensitivity in this context (Supplementary Fig. S7F). Purification of human cells from the secondary PDXs post-therapy and targeted sequencing confirmed their genotype and did not reveal the selection of any further relapse-specific variants (Supplementary Fig. S8A–S8C). The observed differences in drug response could not be accounted for by any consistent changes in L-IC frequency (Supplementary Table S6).

Figure 6.

dRI subclones display decreased sensitivity to commonly used chemotherapeutic drugs. A, Phylogenetic analysis showing clonal relationships in patient 7 based on VAF of leukemic variants shows clear evidence of the isolation of a relapse-fated, dRI clone in dPDX 20. The distances between symbols on the tree were estimated by a nearest neighbour joining method and represent the degrees of relation between them (Minkowski's distance). Circles represent patient samples and triangles represent PDXs; blue represents diagnosis and red represents relapse. dRI-PDX 20 is indicated by a hatched border box. Purified human cells from primary dPDX 7, dRI-PDX 20, and rPDX 5 (representative relapse genetics) xenografts were injected into secondary NSG mice and allowed to engraft. Mice were randomized into 4 groups (with 4 to 5 mice per group) and treated with either saline, dexamethasone, l-asparagine, or vincristine. After 4 weeks of treatment mice were sacrificed and engraftment in the IF, BM, and spleen were analyzed by flow cytometry. Ratio of human chimerism in the BM of drug-treated mice in comparison with saline controls is shown. B, Ratio of human chimerism in the BM of drug-treated mice in comparison with saline controls of dPDX, dRI-PDX, and rPDX of patient 1. C, Representative flow plots of the CD19 and CD33 immunophenotype of dPDX and dRI-PDX dexamethasone-treated mice from patient 1. Lines represent mean and SD. Only significance between dPDX and dRI-PDX are shown. **, P < 0.01; ***, P < 0.001; ****, P < 0.0001; unpaired two-sided t tests.

Figure 6.

dRI subclones display decreased sensitivity to commonly used chemotherapeutic drugs. A, Phylogenetic analysis showing clonal relationships in patient 7 based on VAF of leukemic variants shows clear evidence of the isolation of a relapse-fated, dRI clone in dPDX 20. The distances between symbols on the tree were estimated by a nearest neighbour joining method and represent the degrees of relation between them (Minkowski's distance). Circles represent patient samples and triangles represent PDXs; blue represents diagnosis and red represents relapse. dRI-PDX 20 is indicated by a hatched border box. Purified human cells from primary dPDX 7, dRI-PDX 20, and rPDX 5 (representative relapse genetics) xenografts were injected into secondary NSG mice and allowed to engraft. Mice were randomized into 4 groups (with 4 to 5 mice per group) and treated with either saline, dexamethasone, l-asparagine, or vincristine. After 4 weeks of treatment mice were sacrificed and engraftment in the IF, BM, and spleen were analyzed by flow cytometry. Ratio of human chimerism in the BM of drug-treated mice in comparison with saline controls is shown. B, Ratio of human chimerism in the BM of drug-treated mice in comparison with saline controls of dPDX, dRI-PDX, and rPDX of patient 1. C, Representative flow plots of the CD19 and CD33 immunophenotype of dPDX and dRI-PDX dexamethasone-treated mice from patient 1. Lines represent mean and SD. Only significance between dPDX and dRI-PDX are shown. **, P < 0.01; ***, P < 0.001; ****, P < 0.0001; unpaired two-sided t tests.

Close modal

Unexpectedly, in 2 of 5 patient samples used for drug testing, we found that the dRI subclones exhibited immunophenotypic plasticity as compared with the predominant diagnosis clone when exposed to drugs. In the secondary recipients of dRI-PDX of patient 1 and dRI-PDX of patient 7, both patients harboring DUX4 translocations at diagnosis, there was the emergence of a distinct CD33+CD19dim/− population upon dexamethasone treatment but not in saline controls or upon treatment with the other two drugs (Fig. 6C; Supplementary Fig. S9A and S9B). This population was very rare or not observed in the primary recipients. The CD33+CD19dim/− cells resembled myeloid cells with respect to their cell surface marker expression, size, and granularity (Supplementary Fig. S9C). An immunophenotypic shift toward the myeloid lineage especially with steroid challenge has previously been reported in ERG/DUX4 patients, but our finding of a subclone-specific switch is of interest (34). The emergence of CD33+CD19dim/− cells was not due to the outgrowth of a different subclone, as these cells were genetically identical to their CD19+ counterpart (Supplementary Fig. S9D and S9E). We speculate that this propensity for immunophenotypic plasticity might be linked to chromatin remodeling and the ability to evade therapy, as other studies have shown (35–38). Overall, our data provide strong evidence that prior to exposure to chemotherapy, dRI subclones already possess distinct preexisting functional properties including differential sensitivity to standard chemotherapy agents.

Transcriptional Analysis of dRIs Reveals Metabolic Rewiring and Enrichment of a Stem-Cell State in Progression to Relapse

To gain mechanistic insight into the molecular pathways present in dRIs, RNA-seq analysis was performed on cells from dRI-PDX, dPDX, or rPDX of 4 patients (patients 1, 4, 6, 7; n = 14 dPDXs, 15 dRI-PDXs, and 13 rPDXs; Supplementary Table S7). This analysis confirmed the placement of dRIs as intermediates between diagnosis and relapse, sharing transcriptional programs with both timepoints (Supplementary Fig. S10A). Given the distinct B-ALL subtypes/cytogenetics of the patient samples analyzed, as expected only a few (n = 23) differentially expressed genes reached significance when comparing dPDX versus dRI-PDX across all samples (Supplementary Table S7). Surprisingly, one of the genes upregulated in the dRI-PDX is asparaginase (ASPG), a catalytic enzyme that hydrolyzes asparagine to aspartic acid, albeit nonrobustly in human cells. The human ASPG protein can display cytotoxic activity in human leukemic cell lines (39), suggesting that the dRI-PDX cells may have altered their response to cytotoxic stress or metabolic requirements, thereby explaining why some of the dRI clones are less sensitive to l-asparaginase treatment than dPDX clones.

As pathway enrichment analysis is more sensitive than differential gene expression for finding differences between populations, we undertook this approach to uncover significantly differentially enriched pathways that were shared among all patients. Gene set enrichment analysis (GSEA) comparing the pathways present in dRI-PDX versus dPDX revealed that most pathways significantly enriched in dRI-PDX (FDR q-value ≤ 0.05) were also present or further enriched in the comparison of rPDX to dPDX; these are termed dRI-PDX/rPDX common pathways (Fig. 7A and B; Supplementary Fig. S10B–S10D; Supplementary Table S7). Network analysis revealed that these pathways were centered on a metabolic signature composed of genes involved in cellular and mitochondrial metabolism including amino acid metabolism, tricarboxylic acid cycle, oxidative phosphorylation, mitochondrial translation and transport, and lipid metabolism (Fig. 7A and B; Supplementary Fig. S10B–S10D). In concordance with a previous study on ALL relapse (26), the central regulator of growth and metabolism, mTOR, was enriched in dRI and relapse clones (Supplementary Fig. S10B; Supplementary Table S7). Pathways identified as uniquely enriched in rPDX versus dPDX (rPDX-unique) included a large network of highly interconnected nodes involved in cell-cycle regulation such as cell-cycle checkpoints, DNA replication, DNA repair, and microtubule organization (Fig. 7A; Supplementary Table S7).

Figure 7.

dRI subclones share a common metabolic and stem cell profile. A, Plot showing the normalized enrichment score (NES) for the top differentially enriched gene sets (FDR q value ≤ 0.05) of dRI-PDX unique, dRI-PDX/rPDX common, and rPDX unique groups from GSEA of the following comparisons: dPDX vs. dRI-PDX and dPDX vs. rPDX. B, Heat maps showing the expression of leading-edge genes (subset of genes found in the ranking at or just before the maximal enrichment score in GSEA) for selected gene sets enriched in dRI-PDX and rPDX from enrichment map in A. Relative expression was generated from variance stabilized normalized counts C, dPDX, dRI-PDX, and rPDX from Patient 1, Patient 4, and Patient 7 were stained with MitoTracker and CellROX dyes and analyzed by flow cytometry. MFI for each sample and dye is represented as ratio to dPDX for each patient (Patient 1: dPDX n = 5, dRI-PDX n = 5, rPDX = 4; Patient 4: dPDX n = 5, dRI-PDX n = 4, rPDX = 4; Patient 7: dPDX n = 5, dRI-PDX n = 5, rPDX = 5). D, Immunostaining analysis for TOMM20, MRPS18B, and SOD1 in dPDX, dRI-PDX, and rPDX cells from Patients 1, 3, and 7. The Integrated Density (IntDen = Area × MFI) for 40 to 50 cells from each clone was analyzed using Fiji. The mean for each clone was normalized and calculated as a ratio to the dPDX for each patient separately. Representative images for Patient 7 are shown. (Patient 1: dPDX n = 1, dRI-PDX n = 2, rPDX = 1; Patient 7: dPDX n = 2, dRI-PDX n = 2, rPDX = 1; Patient 3: dPDX n = 2, dRI-PDX n = 2, rPDX = 1). E, GSEA enrichment plots from the following comparisons: dPDX vs. dRI-PDX (n = 4 pts); dPDX vs. rPDX (n = 4 pts); and diagnosis (Dx) versus relapse (Rel) patient samples from our cohort were generated for mitochondrial translation and oxidative phosphorylation gene sets. F, Bar plot of the aggregated GSVA scores for B-cell genes and HSC genes in each sample. GSVA scores for samples in each category were summed and scaled from 0 to 1. The numbers above the bars represent how many times the observed score was higher than random scores obtained in 1,000 permutations using a list of 1,000 random genes. G, Schematic diagram of dRI with altered metabolic and stem-cell programs preexisting in diagnosis patient samples that survive chemotherapy and seed the relapse disease. *, P < 0.05; **, P < 0.01; ***, P < 0.001; unpaired two-sided t test.

Figure 7.

dRI subclones share a common metabolic and stem cell profile. A, Plot showing the normalized enrichment score (NES) for the top differentially enriched gene sets (FDR q value ≤ 0.05) of dRI-PDX unique, dRI-PDX/rPDX common, and rPDX unique groups from GSEA of the following comparisons: dPDX vs. dRI-PDX and dPDX vs. rPDX. B, Heat maps showing the expression of leading-edge genes (subset of genes found in the ranking at or just before the maximal enrichment score in GSEA) for selected gene sets enriched in dRI-PDX and rPDX from enrichment map in A. Relative expression was generated from variance stabilized normalized counts C, dPDX, dRI-PDX, and rPDX from Patient 1, Patient 4, and Patient 7 were stained with MitoTracker and CellROX dyes and analyzed by flow cytometry. MFI for each sample and dye is represented as ratio to dPDX for each patient (Patient 1: dPDX n = 5, dRI-PDX n = 5, rPDX = 4; Patient 4: dPDX n = 5, dRI-PDX n = 4, rPDX = 4; Patient 7: dPDX n = 5, dRI-PDX n = 5, rPDX = 5). D, Immunostaining analysis for TOMM20, MRPS18B, and SOD1 in dPDX, dRI-PDX, and rPDX cells from Patients 1, 3, and 7. The Integrated Density (IntDen = Area × MFI) for 40 to 50 cells from each clone was analyzed using Fiji. The mean for each clone was normalized and calculated as a ratio to the dPDX for each patient separately. Representative images for Patient 7 are shown. (Patient 1: dPDX n = 1, dRI-PDX n = 2, rPDX = 1; Patient 7: dPDX n = 2, dRI-PDX n = 2, rPDX = 1; Patient 3: dPDX n = 2, dRI-PDX n = 2, rPDX = 1). E, GSEA enrichment plots from the following comparisons: dPDX vs. dRI-PDX (n = 4 pts); dPDX vs. rPDX (n = 4 pts); and diagnosis (Dx) versus relapse (Rel) patient samples from our cohort were generated for mitochondrial translation and oxidative phosphorylation gene sets. F, Bar plot of the aggregated GSVA scores for B-cell genes and HSC genes in each sample. GSVA scores for samples in each category were summed and scaled from 0 to 1. The numbers above the bars represent how many times the observed score was higher than random scores obtained in 1,000 permutations using a list of 1,000 random genes. G, Schematic diagram of dRI with altered metabolic and stem-cell programs preexisting in diagnosis patient samples that survive chemotherapy and seed the relapse disease. *, P < 0.05; **, P < 0.01; ***, P < 0.001; unpaired two-sided t test.

Close modal

To validate the metabolic signature observed in dRIs that was often further enhanced at relapse, we performed flow cytometry and immunostaining analysis in the PDX from five patient samples (Fig. 7C and D; Supplementary Fig. S10E and S10F). This analysis confirmed the convergence of similarities in metabolic rewiring of dRI-PDX/rPDX in comparison with dPDX in all 5 patients with B-ALL. This functional validation showed that there was an increase in total mitochondrial mass, with individual mitochondria having similar mitochondrial membrane potential in most of the patients analyzed. However, despite this increase in absolute mitochondria abundance, the levels of reactive oxygen species (ROS) were unexpectedly found to be lower in dRI-PDX and rPDX cells as compared with dPDX cells. The reduced ROS levels are suggestive of the presence of increased antioxidant defense upon progression toward relapse. This interpretation was concordant with our gene-expression data that showed an enrichment of ROS-defense and peroxisomal activity genes (i.e., catalase), which degrade toxic hydrogen peroxide and metabolize drugs, in the dRI-PDX and rPDX samples (Fig. 7A and C; Supplementary Fig. S10B, S10D–S10F; Supplementary Table S7). Significant enrichment for chromatin remodeling and cell stress response (such as the UPR) were also identified as dRI-PDX/rPDX common pathways (Fig. 7A; Supplementary Fig. S10B–S10D). The chromatin remodeling pathways including expression of histone variants and isoforms were likely instrumental in the plasticity observed in the progression to relapse, including the immunophenotypic plasticity we observed in dRI-PDXs of patients 1, 7, and 9 (Figs. 4A and 5C), whereas the expression of the cell stress response pathways identified could contribute to enhanced survival of dRI subclones. To independently validate the results obtained from the PDX-identified pathways, we directly evaluated bulk diagnosis and relapse patient samples. We found enrichment (FDR q-value ≤ 0.05) of pathways involved in metabolism, mitochondrial regulation, and cell cycle at relapse (Fig. 7E; Supplementary Fig. S10G; Supplementary Table S7). Furthermore, we found significant enrichment of a signature present at minimal residual disease (18) in dRI-PDX and rPDX compared with dPDX, lending further support for the relevance of the dRI-PDX/rPDX pathways for understanding the mechanisms of relapse disease initiation and their presence prior to chemotherapeutic challenge (Supplementary Fig. S10H). Collectively, our transcriptional analysis has for the first time uncovered the utilization of chromatin remodeling, stress responses, and metabolic pathways in dRI subclones whose activity could serve to protect them during chemotherapy treatment and contribute to their ability to further progress to relapse disease. As such, these newly identified pathways represent rich areas to investigate for new therapeutic strategies to target dRI specifically.

To investigate whether the broader cellular state of the dRIs contributes to their functional differences, we performed gene set variation analysis (GSVA) comparing the transcriptomic profiles of representative clones with normal hematopoietic cell populations isolated from human umbilical cord blood using two independent datasets, one from our own newly generated data and the other published previously (ref. 40; Supplementary Fig. S11A; Supplementary Table S7). dRI-PDX and rPDX exhibited a transcriptome profile significantly enriched for hematopoietic stem cell (HSC) genes and a slight reduction in B-cell genes as compared with dPDX (Fig. 7F; Supplementary Fig. S11B and S11C). Depletion of B-cell genes and enrichment of HSC genes at relapse were also observed in the bulk patient samples from our cohort (Fig. 7F). Furthermore, the leading-edge genes of common pathways enriched in both dRI-PDX and rPDX are upregulated in HSCs as compared with lymphoid cells (Supplementary Fig. S11D). The enrichment of mitochondrial metabolism and stemness at relapse were also validated in the larger Waanders and colleagues cohort of paired diagnosis-relapse patient B-ALL samples (7) encompassing several genetic subtypes (Supplementary Fig. S11E and S11F). These findings are in line with previous studies in the normal and leukemia stem-cell fields (25, 41, 42) where several metabolic and stress response pathways such as those identified here (mitochondria metabolism, UPR, antioxidant defense) have been described to be crucial to the maintenance of stem-cell homeostasis and function. In addition, enrichment of stemness signatures is also a hallmark of high-risk B-ALL (43). Thus, our findings report a link for the acquisition of HSC stemness properties in combination with metabolic rewiring in dRI as part of progression to relapse in lymphoid B-ALL (Fig. 7G).

Our study provides new insight into the leukemogenic process of human B-ALL through a deep analysis of the functional and molecular properties of genetically diverse diagnosis subclones isolated through xenografting from pediatric and adult patients with B-ALL. By combining xenografting with sequencing, broad clonal structures were unambiguously demonstrated and additional subclones were uncovered. We identify relapse-fated dRI subclones, prior to chemotherapeutic exposure, with the capacity for clonal propagation and leukemia initiation in B-ALL that are both genetically and transcriptionally related to the relapse. Our data extends prior B-ALL studies that relied on genetic analysis and computational methods of bulk leukemia samples without functional studies, and that could only infer the presence and survival of minor subclones that were ancestral to the relapse driving clone and could not speak to the functional properties of these cells at diagnosis (4, 5). As shown by Waanders and colleagues (7), although our xenografting approach substantiates the predicted lineage relationships of the traditional genetic analysis of bulk tissues, it also extends our understanding of leukemia evolution by revealing a richer diversity of branching evolution and additional subclones. Our functional studies of subclones analyzed either alone or in the context of competitive repopulation demonstrate that dRI subclones can be found to display increased clonal repopulation kinetics, immunophenotypic plasticity, and differences in leukemic dissemination such as dissemination to the CNS, a known disease sanctuary site in patients. Importantly, dRI showed an intrinsic difference in tolerance to standard induction chemotherapeutic agents. Moreover, dRI subclones were typically very minor latent subclones in the diagnostic sample, suggesting that they are restrained within the patient in comparison to the predominant subclone, perhaps due to dormancy, timing, subclonal cross-talk, or slow cycling, although this was not directly measured in our studies. Our study resolved the theoretical question as to whether functional properties of relapse clones are present prior to chemotherapeutic challenge or are induced by chemotherapy. Although some mutations not present at diagnosis are recurrently acquired after chemotherapy (e.g., USH2A, NT5C2), our study documents the preexistence of dRI subclones with intrinsic tolerance to clinical drugs. This data provides direct evidence that the Luria–Delbrück principle that resistance in a cell population may be intrinsic occurs in human leukemia (44). We speculate that the dRI subclones arise stochastically and these partially evolved subclones then become selected for by chemotherapy. Gene-expression and pathway analysis showed that even at this early stage where relapse fate arises, there was considerable metabolic rewiring within dRI subclones that was often further enhanced in relapse clones. Recurrent pathways were identified across a spectrum of different patients with B-ALL with distinct genetic drivers, suggestive of convergent evolution. Broadly, we found that adult and pediatric dRI had similar properties in terms of their functional xenografting properties, therapy response, transcriptional pathways, and metabolic rewiring, also supporting the concept of convergent evolution. Pathways regulating mitochondrial dynamics and proteostasis are particularly relevant because they are linked to survival under stress, processes that may underpin the ability of dRI cells to adapt and survive chemotherapy. Further evolution and mutagenesis of such subclones during the treatment and/or remission phases ultimately results in disease recurrence. Expression of the dRI-PDX/rPDX common pathways converged on an HSC signature, shared by both dRI and relapse cells, analogous to findings from recent studies we have undertaken from purified AML stem-cell populations (24, 45). The finding of stemness signatures already present in dRI subclones, together with the identified metabolic changes, provides the molecular basis to explain why dRI subclones may both survive therapy and possess the regenerative capacity to initiate disease relapse after a period of dormancy. Future studies that explore the transcriptome of dRI subclones from across an even wider spectrum of B-ALL patient samples would give insight into how the timing of acquisition of these molecular pathways leads to relapse.

The ability to isolate and characterize dRI subclones provides an important first step in understanding the basis for therapy resistance and clonal propagation. Such subclones isolated in this way are not simply an in silico depiction; rather, they can be viably preserved cells or serially propagated xenografts for future studies. Whole-genome sequencing, methylation, and chromatin accessibility studies could be undertaken to build upon our transcriptomic analysis, thereby exploring more deeply the mechanisms that drive their functional properties. The signatures developed from isolated dRI clones could yield biomarkers to improve the classification of patients who are at increased risk of relapse and to better monitor residual disease. The immunophenotypic plasticity we linked to different genetic subclones points to the need to understand the breadth of cell-surface phenotypes present to ensure that all leukemic subclones are properly tracked during flow cytometry–based residual disease monitoring. Finally, further investigation of the dRI transcriptomic profile and metabolic rewiring may be used to uncover the vulnerabilities of dRI subclones, resulting in new therapeutic targets. Improved eradication of dRIs during early treatment phases before the subclones evolve would prevent progression to more aggressive, therapy-resistant disease. Prior studies from us and others have shown that xenografting of a wide spectrum of primary patient samples provides a powerful tool to evaluate novel therapeutics and develop biomarkers (23, 25, 46, 47). Our study suggests that extending the xenograft-based drug development paradigm by including genetic analysis to uncover subclonal responses to drug treatment will open up avenues to evaluate whether relapse-fated clones are effectively targeted.

Patient Samples

Patient samples were obtained at diagnosis and relapse from 6 adult patients and 8 pediatric patients with B-ALL according to preestablished guidelines approved by the Research Ethics Board of the University Health Network and the St. Jude Institutional Review Board, respectively, and were conducted in accordance with recognized ethical guidelines. Adult samples were collected at the Princess Margaret Cancer Centre (Toronto, Ontario, Canada) and pediatric samples were collected at St. Jude Children's Research Hospital (Memphis, TN). Written informed consent was obtained from all patients or patient families. All samples were frozen viably and stored long term at −150°C. Samples were selected retrospectively based on sample and paired sample availability.

PDX Generation and Limiting Dilution Analysis

Twenty-nine clinical samples obtained from the 14 patients (1 patient having two relapse samples) were stained with the following antibodies: anti-CD19 PE (BD Biosciences, clone 4G7), anti-CD3 FITC (BD Biosciences, clone SK7) or anti-CD3 APC (Beckman Coulter clone UCHT11), anti-CD45 APC (BD Biosciences, clone 2D1) or anti-CD45 FITC (BD Biosciences, clone 2D1), and anti-CD34 APC-Cy7 (clone 581). Each sample was sorted on a FACSAria III (BD Biosciences) for leukemic blasts (CD19+CD45dim/−) and T cells (CD3+ CD45hi). NSG mice were bred according to protocols established and approved by the Animal Care Committee at the University Health Network. Eight- to 12-week-old mice were sublethally irradiated at 225 cGy 24 hours before transplant. Only female mice were used in these studies. Intrafemoral injections of 10 to 250,000 sorted leukemic blasts were performed as described previously (48). Mice were sacrificed 20 to 30 weeks post-transplant or at the onset of disease symptoms. Human cell engraftment in the injected femur, bone marrow (noninjected bones, left tibia, right tibia, left femur), spleen, and CNS were accessed using human specific antibodies for CD45 (PE-Cy7, BD Biosciences, clone HI30; v500, BD Biosciences, clone HI30), CD44 (PE, BD Biosciences, clone 515; FITC, BD Biosciences, clone L178), CD3 (APC, BD Biosciences, clone UCHT1), CD19 (PE-Cy5, Beckman Coulter, clone J3-119), CD33 (PE-Cy7, BD Biosciences, clone P67-6; APC, BD Biosciences, P67-6), and CD34 (APC-Cy7, BD Biosciences, clone 581) analyzed on an LSRII (BD Biosciences). Mice were considered to be engrafted when >0.1% of cells in the injected femur were positive for one or more human B-ALL specific cell surface markers (CD45, CD44, CD19, CD34). The confidence intervals for the frequency of L-ICs was calculated using ELDA software (49).

WES

DNA from the adult samples was extracted from sorted leukemic blasts (CD19+CD45dim/−, purity >90%) and sorted T cells (CD3+CD45hi, purity >90%) using the Gentra Puregene Cell Kit (Qiagen). T cells from patients 2, 3, and 5 were whole-genome amplified (REPLI-g Mini Kit, Qiagen) due to limited material. DNA from the pediatric samples with >90% leukemic blasts was extracted from diagnosis, remission, and relapse samples using phenol–chloroform. Exomes were captured using the TruSeq Exome Library Prep Kit (67 Mb, 1 μg DNA input) or the Nextera Rapid Capture Expanded Exome (62 Mb, 50 ng DNA input; Illumina). Paired-end sequencing was performed with the HiSeq 2500 genome sequencer (Illumina). The data was mapped to human reference genome hg19 and variant calling was performed using the Bambino variant detector as described previously (50). Briefly, leukemia and germline files were combined by the program and aligned against the reference genome. Putative sequence variants including SNVs and indels were detected by running the variation detection module of Bambino. The output contained detailed read counts for each variant with columns for tumor/normal status, allele, and strand. Variants were not filtered for coverage prior to combination with targeted-sequencing results.

Copy-Number Analysis

Patient copy-number aberrations were determined using SNP6.0 microarrays according to manufacturer's instructions (Affymetrix). Data was analyzed as described previously (51) using optimal reference normalization (52) and circular binary (53, 54) segmentation with Genotyping Console (Affymetrix) and dCHIP (build Apr 2010; ref. 55). Detection of loss of heterozygosity and allelic ratios were performed using Nexus 7.5.2 software (BioDiscovery Inc). All segments were manually curated.

Recovery of Human Cells and DNA Isolation from Xenografts

Cells from the injected femur (IF), BM, and spleen were frozen viably after sacrifice. The IF and BM of mice engrafted with >10% human cells were combined. These cells as well as cells from diagnosed PDX spleens were depleted of mouse cells using the Miltenyi Mouse Cell Depletion Kit (Miltenyi Biotec; samples with >20% engraftment) or by cell sorting with human CD45 and human CD19 and or CD34 cell-surface antibodies to a purity of >90% as determined by postprocessing flow cytometry. CNS cells from mice with greater than 60% engraftment were used directly for DNA isolations. DNA was isolated using the QIAamp DNA Blood Mini or Micro Kit (Qiagen).

Targeted Sequencing

All somatic SNVs and indels identified by WES were validated in the patient samples using NimbleGen SeqCap Target Enrichment according to the manufacturer's instructions (Roche, NimbleGen). Library preparation was completed using 250 to 500 ng of DNA using the NEXTflex DNA-SEQ Library Prep Kit (BiooScientific) with NEXTflex-96 DNA Barcodes (BiooScientific). Sequencing was performed on a HiSeq 2500 genome sequencer to a mean coverage >350× for patient samples and >200× for PDXs.

Targeted-Sequencing Data Analysis

Final patient variant calls used the combined results of WES and capture validation. Variants were filtered out if the VAF in the germline was greater than 10% or if there was a dbSNP frequency of greater than 1%. Variants were classified based on the VAFs in the bulk patient diagnosis and relapse samples as: preserved variants (VAF > 30% in both diagnosis and relapse samples, or preserved between samples); diagnosis-specific variants [present at diagnosis (>1%) and absent at relapse (<1%)]; latent variants (present at diagnosis with VAF < 30% and increasing at relapse); relapse-specific variants [absent at diagnosis (<1%) and present at relapse (>1%)]. For xenograft analysis, variants with less than 5× coverage or uncovered in numerous xenografts were removed. Results were analyzed by visualization in heat maps (i.e., Fig. 1C). Phylogenetic analysis showing genetic relationships of patient samples and xenografts were estimated using Minkowski's distance calculated from the VAFs and represented visually using a nearest neighbor joining algorithm. Genetic concordance between different tissues of the same xenograft were determined by visual assessment by three independent and blinded individuals. Discordance was called only when all three investigators were in agreement.

Generation of Mutational Trees from Patient Samples and Xenografts

Two independent computational analyses were performed—first for patient-only tissue samples, and then for patient samples augmented with xenografts (BM and spleen)—using an early version of the PairTree algorithm (Wintersinger and colleagues, in preparation). PairTree uses variant read counts to estimate the posterior probability distribution over four possible ancestral relationships between every variant pair (A, B). The four ancestral relationships are as follows: variant A and variant B occurred in the same cells, such that no cell possessed one variant but not the other; variant A is ancestral to variant B, such that some cells have A but not B; variant B is ancestral to variant A, such that some cells have B but not A; or neither is the ancestor of the other, such that A and B are on different branches of the evolutionary tree.

To permit temporal ordering of mutations, the infinite sites assumption was made, such that variant A could never be the ancestor of variant B if A's cancer cell frequency (CCF) was lower than B's in any sample, after their CCFs were estimated from each mutation's VAF. To simplify the estimation of CCF from VAF, variants were discarded if they lay in CNA-affected regions determined by SNP6.0 analysis of the patient samples, or if their VAFs were suggestive of an uncalled CNA region. X-chromosome variants in male patients had their VAFs halved when estimating CCF to correct for their haploid nature. Distinct mutational populations were defined by semiautomatically clustering variants, with variants clustered together when their mutual “A and B occurred in the same cell” probabilities were high.

For each cluster, a representative “supervariant” was created. To compute a precise estimate of supervariant CCF, we summed the read counts of all variants within the supervariant's cluster, then computed the pairwise ancestral relationship probabilities between supervariants. These probabilities guided expert-driven tree construction, in which trees were built according to how well their implied pairwise relationships between clusters matched the computed pairwise probabilities. Confirmation that these expert-driven trees yielded data likelihoods as good as or better than fully automated tree reconstructions produced by a more mature version of the PairTree algorithm (Wintersinger and colleagues, in preparation) was performed. The tree-constrained lineage frequencies of each population in each sample were then computed using PairTree. The variants included in the populations defined by the patient-only analysis were then compared to the variants included in the populations defined for the patient and xenograft combined analysis to reveal the overlap, differences, and additional clonality described with inclusion of the xenograft samples.

RNA-seq

For PDX

Cell pellets for PDX RNA extraction were frozen in 1 mL of TRIzol (Invitrogen) and kept at −80°C. Total RNA was purified by phenol/chloroform and integrity and concentration were verified on a Bioanalyzer Pico Chip (Agilent Technologies). cDNA conversion was performed using SMART-Seq v4 Ultra Low Input RNA Kit for Sequencing (Takara; 1 ng total RNA input) and libraries were prepared using Nextera XT DNA Library Preparation Kit (Illumina; 1 ng input of cDNA). Equimolar quantities of libraries were pooled and sequenced 4 cDNA libraries per lane on a High Throughput Run Mode Flowcell with v4 sequencing chemistry on the Illumina HiSeq 2500 following manufacturer's protocol generating paired-end reads of 126-bp in length to reach depth of 65 million reads per sample.

For Patient Samples

RNA was extracted from sorted leukemic blasts (CD19+CD45dim/−, purity >90%) using RNeasy Micro Kit (Qiagen) for adult patient samples. From these samples approximately 10 ng of total RNA was processed using the SMART cDNA synthesis protocol including SMARTScribe Reverse Transcriptase (Clontech) as per the manufacturer's instructions. The amplified cDNA was subject to automated Illumina paired-end library construction using the NEBNext paired-end DNA sample Prep Kit (NEB) and libraries were sequenced on HiSeq2000 (Illumina) to an average of approximately 161 million Chastity-passed paired reads of 75 bp in length per sample. For pediatric samples RNA was extracted from patient samples with >90% leukemic blasts using TRIzol (Life Technologies) and quality and quantity assessed by Qubit (Thermo Fisher Scientific) and RNA6000 Chip (Agilent). One microgram of RNA was used for library preparation with TruSeq RNA Library Prep Kit v2 (Illumina) and 2 × 100 bp paired-end sequencing was performed on a HiSeq 2500 (Illumina). Patient RNA-seq samples were aligned against the human genome (hg19) using STAR v2.3 with default parameters. All PDX samples were aligned with STAR v2.5.2b (56) against the human genome build version GRCh38 and the ensembl v90 gtf file. Default parameters were used except chimeric segments were screened with a minimum size 12, junction overlap 12, and segment reads gap maximum 3; splice junction overlap 10, aligned mates gap maximum 100,000, aligned intron maximum 100,000 and alignSJstitchMismatchNmax 5 -1 5 5. For both patient samples and PDX, transcript counts were obtained using HTSeq v0.7.2 (57). Data was library size normalized using RLE, followed by a variance stabilizing transformation using DESeq2 v1.22.1 (58). Principal component analysis plots were generated on a per sample basis using the top 1,000 variable genes. For downstream visualization, differential expression, and pathway analysis, the mean expression of each sample clone condition was utilized. On a per patient level, differentially expressed genes were identified. Genes with adjusted P value of < 0.05 and absolute log2 fold change of > 1 were considered significant. All visualizations were generated using R v3.5.1 and the pheatmap v1.0.10 and lattice v0.20–38 packages and ggplot2 v3.1.0 packages.

Pathway Enrichment Analysis and Visualization

Pathway enrichment analysis and visualization was performed as described previously (59). Briefly, a score to rank genes from top upregulated to downregulated was calculated using the formula −sign(logFC) × −log10(P). The rank file from each comparison was used in GSEA (http://software.broadinstitute.org/gsea/index.jsp) using 2,000 permutations and default parameters against a pathway database containing Msigdb c2 and c3, NCI, IOB, NetPath, HumanCyc, GO BP, Reactome, and Panther (http://baderlab.org/GeneSets, version June 2018). GSEA progressively examines genes from the top to the bottom of the ranked list, increasing the enrichment score (ES) if a gene is part of the pathway and decreasing the score otherwise. These running sum values are weighted, so that enrichment in the very top-ranking (and bottom-ranking) genes is amplified, whereas enrichment in genes with more moderate ranks is not amplified. The ES is calculated as the maximum value of the running sum and normalized relative to pathway size, resulting in a normalized ES (NES) that reflects the enrichment of the pathway in the list. Positive and negative NES values represent enrichment at the top and bottom of the list, respectively. A permutation-based P value is computed and corrected for multiple testing to produce a permutation-based FDR q-value that ranges from zero (highly significant) to one (not significant). EnrichmentMap version 3.1.0 in Cytoscape 3.7.0 was used to visualize enriched gene sets at FDR q-value ≤ 0.05 for each comparison with a Jaccard coefficient set to 0.375. The Enrichment Map software takes as input a text file containing pathway enrichment analysis results (from GSEA) and another text file containing the pathway gene sets used in the original enrichment analysis (http://baderlab.org/GeneSets, version June 2018). An enrichment map is a network representing overlaps among enriched pathways with pathways represented as circles (nodes) that are connected with lines (edges) to other pathways with overlapping genes. The network layout and clustering algorithm Autoannotate app in Cytoscape was used to automatically display and group similar pathways as major biological themes. Pathways differentially enriched in dPDX versus dRI-PDX and/or dPDX versus rPDX clones were classified as: dRI-PDX-unique [pathways significantly enriched (FDR q-value ≤ 0.05) in dRI-PDX in comparison with dPDX], dRI-PDX/rPDX common [pathways significantly enriched (FDR q-value ≤ 0.05) in both dRI-PDX and rPDX in comparison with dPDX] or rPDX-unique [pathways significantly enriched (FDR q-value ≤ 0.05) in rPDX in comparison with dPDX].

Isolation of HSCs and B Cells from Human Cord Blood Samples

Human cord blood samples were obtained with informed consent from Trillium and Credit Valley Hospital according to procedures approved by the University Health Network Research Ethics Board. Mononuclear cells were obtained by centrifugation on Lymphoprep medium (StemCell Technologies) followed by red blood cell lysis using ammonium chloride (StemCell Technologies). Human CD34+ and CD34 CB cells were separated using CD34 Microbeads Kit (Miltenyi Biotec) according to manufacturer's protocol and stored at −150°C. Cells were stained with antibodies (all from BD Biosciences, unless stated otherwise): FITC–anti-CD45RA (1:50, 555488), PE–anti-CD90 (1:50, 555596), PECy5–anti-CD49f (1:50, 551129), V450–anti-CD7 (1:33.3, 642916), PECy7–anti-CD38 (1:100, 335790), APC–anti-CD10 (1:50, 340923), APCCy7–anti-CD34 (1:200, custom made by BD Biosciences), and sorted on FACS Aria III (Becton Dickinson), consistently yielding >95% purity. HSCs were sorted on the basis of the following markers: CD34+CD38CD45RACD90+ CD49f+ from CD34+ CB cells as described previously (60). B cells were sorted as CD34CD38+CD19+CD33CD3CD56.

GSVA

Read counts from 5 B-cell and 3 HSC-sorted populations from normal human umbilical cord blood were normalized using the trimmed of M values using the edgeR_3.24.3 R package, and gene differential expression was calculated using the quasi likelihood ratio method and included cord blood batch correction in the experimental design. The top 1,000 genes enriched in B cells (B-cell genes) and the top 1,000 genes enriched in HSCs (HSC genes) were selected to be used as the reference signature. Two mixed profiles were used: dPDX, dRI-PDX and rPDX clones and 14 bulk diagnosis and relapse patient samples. The mixed profiles read counts were variance stabilized and library size normalized using DESeq_1.34.1. The gsva() function from the R GSVA_1.30.0 package was employed using a Gaussian kernel to estimate enrichment of reference signatures in each sample of the mixed profile dataset. GSVA score for each patient and each category were plotted on a box plot and a strip chart. GSVA scores were summed for each mixed profile category and standardized from 0 to 1. One thousand permutations with a random gene list of size 1,000 were performed on the mixed profile and percentages calculated to indicate how many times the observed score was higher than the random scores. These results were confirmed using the HSC and B-cell expression profiles from Novershtern and colleagues (ref. 40; Supplementary Table S7).

Transcriptomic Validation Experiments

Staining for mitochondria content, ROS, mitochondrial membrane potential, and mitochondrial superoxide levels was performed by incubating thawed PDX cells at 37°C with 1 mmol/L of MitoTrackerGreen (M7514; dilution 1:10,000); 5 μmol/L CellROX deep red (C10422; dilution 1:500); 1 μmol/L TMRE (T668; dilution 1:1,000), or 5 mmol/L MitoSOX Red (M36008; dilution 1:1,000) following the manufacturer's instructions (Thermo Fisher Scientific) and directly analyzed on a BD LSRII cytometer. Mean fluorescence intensity (MFI) for each sample and dye is represented as ratio to dPDX for each patient. Immunostaining analysis for TOMM20, MRPS18B, and SOD1 were performed on PDX cells. Briefly, cells were spun onto poly-l-lysine (Sigma)–coated slides (Ibidi, 200 × g, 10 minutes), fixed with 4% paraformaldehyde (Sigma) and permeabilized with 0.5% Triton (Sigma) before blocking (PBS, 10% FBS, 5% BSA). Slides were incubated with primary antibodies (TOMM20: ab56783, dilution 1/200; SOD1: ab8866, dilution 1/100; and MRPS18B: ab191891, dilution 1/200) in blocking solution overnight at 4°C. Secondary anti-mouse AF568, anti-rabbit AF488, and anti-sheep AF647 (Invitrogen, 1:400) antibodies were added (PBS, 0.025% Tween, Sigma, 1.5 hours, room temperature). After washing, nuclei were stained with 1 μg/mL DAPI (Invitrogen) and slides were mounted (Fluoromount G, Invitrogen). Images were captured by a Zeiss LSM700 Confocal (oil, 63×/1.4NA, Zen 2012) and analyzed with ImageJ/Fiji. The Integrated Density (IntDen = Area × MFI) for 40 to 50 cells from each clone was analyzed. The mean for each clone was normalized and calculated as a ratio to the dPDX for each patient separately.

Secondary Transplantations for Drug Assays and Limiting Dilution Assays

Human purified cells from the primary recipients were thawed and transplanted into 8- to 12-week-old female NSG mice sublethally irradiated as described for the primary recipients. The number of mice used for secondary transplantation experiments/drugs was determined by cell and mouse availability and feasibility. Intrafemoral injections of 26,000 to 100,000 leukemic blasts were performed for drug assays and a range of 10,000 to 100,000 leukemic blasts were injected for secondary limiting dilution assays. After 4 weeks, mice were bled from the saphenous vein and human chimerism was evaluated by flow cytometry. Once human engraftment in the peripheral blood reached between 1% and 10%, or after 10 to 14 weeks for those samples in which leukemic blasts did not mobilize to the peripheral blood, mice were randomized and single-agent treatments were started. Dexamethasone (15 mg/kg), l-asparaginase (1,000 kU/kg), and saline were given daily by intraperitoneal injection 5 days a week. Vincristine (0.5 mg/kg) was given once a week by intraperitoneal injection. All four arms of the drug treatment were performed for 4 weeks and mice were sacrificed the following day (vincristine/dexamethasone/saline) or 1 week after the last treatment (vincristine). Analysis of the secondary limiting dilution assay was performed 16 weeks posttransplant. Human cell engraftment in the injected femur, bone marrow, and spleen were assessed using human-specific antibodies for CD45 (FITC, BD Biosciences, clone 2D1; v500, BD Biosciences, clone HI30), CD19 (BD Biosciences, PE, clone 4G7), CD33 (PE-Cy7, BD Biosciences, clone P67-6), CD3 (APC, BD Biosciences, clone UCHT1), CD10 (Qdot605, BD Biosciences, clone HI10A), CD38 (BV421, BioLegend, clone HIT2), and CD34 (APC-Cy7, BD Biosciences, clone 581). Mice were considered to be engrafted when >0.1% of cells in the injected femur were positive for one or more human B-ALL–specific cell-surface markers (CD45, CD44, CD19, CD34). Response to treatment was analyzed as a ratio of human engraftment of drug-treated versus saline-treated mice to eliminate interclonal differences in engraftment levels. Ratio of human engraftment in each individual drug-treated PDX to the average engraftment level of all saline controls was calculated. Lineage stains were performed on xenografts expressing CD33 (APC, BD Biosciences clone P67-6) and CD19 (PE, BD Biosciences, clone 4G7) including CD14 (PE-Cy7, Beckman Coulter, clone 52), CD15 (v450, BD Biosciences, clone MMA), CD10 (Qdot605, BioLegend, clone HIT2), and CD34 (APC-Cy7, BD Biosciences, clone 581). The confidence intervals for the frequency of L-ICs was calculated using ELDA software (49). Statistical analysis was performed using PRISM 6 (GraphPad Software).

FACS from Xenografts

Leukemic cells from primary xenografts were sorted for immunophenotypic populations on a FACSARIAIII (BD Biosciences). Cells from the IF and BM were pooled and stained with CD19, CD34, CD45, CD10, and CD33 and collected at a sort purity of >99%.

Data Availability Statement

The datasets generated during this study are included in this published article and its supplementary information files.

Code Availability Statement

Code used in this study is available at https://www.github.com/morrislab/pairtree

J.S. Danska reports receiving commercial research support from Trillium Therapeutics, Inc. and has ownership interest (including patents) in the same. M.D. Minden reports receiving commercial research support from Kura and has received speakers bureau honoraria from Astellas. C.G. Mullighan is a scientific advisory board member for Illumina, reports receiving commercial research grants from AbbVie, Pfizer, and Loxo Oncology, and has received speakers bureau honoraria from Illumina, Amgen, Pfizer, and Aptitude Health.

J.E. Dick served on the SAB at Trillium Therapeutics, reports receiving a commercial research grant from Celgene, and has ownership interest (including patents) in Trillium Therapeutics Inc. No potential conflicts of interest were disclosed by the other authors.

Conception and design: S.M. Dobson, C.G. Mullighan, J.E. Dick

Development of methodology: S.M. Dobson, J. Wintersinger, I. Grandal, G. Bader, J. Easton, J.S. Danska, Q. Morris, C.G. Mullighan, J.E. Dick

Acquisition of data (provided animals, acquired and managed patients, provided facilities, etc.): S.M. Dobson, L. García-Prat, R.J. Vanner, E. Waanders, J. McLeod, O.I. Gan, I. Grandal, S.Z. Xie, M. Hosseini, S.R. Olsen, G. Neale, S.M. Chan, J. Easton, C.J. Guidos, J.S. Danska, M.D. Minden, C.G. Mullighan, J.E. Dick

Analysis and interpretation of data (e.g., statistical analysis, biostatistics, computational analysis): S.M. Dobson, L. García-Prat, R.J. Vanner, J. Wintersinger, E. Waanders, Z. Gu, M.N. Edmonson, X. Ma, Y. Fan, V. Voisin, M. Chan-Seng-Yue, S. Abelson, M. Rusch, S.M. Chan, G. Bader, J. Zhang, M.D. Minden, Q. Morris, C.G. Mullighan, J.E. Dick

Writing, review, and/or revision of the manuscript: S.M. Dobson, L. García-Prat, R.J. Vanner, J. Wintersinger, E. Waanders, O.I. Gan, S.Z. Xie, S. Abelson, J. Easton, J.S. Danska, M.D. Minden, Q. Morris, C.G. Mullighan, J.E. Dick

Administrative, technical, or material support (i.e., reporting or organizing data, constructing databases): S.M. Dobson, L. García-Prat, J. Wintersinger, E. Waanders, J. McLeod, I. Grandal, D. Payne-Turner, P. Gupta, Y. Shao

Study supervision: Q. Morris, C.G. Mullighan, J.E. Dick

We thank the patients and their families and physicians who made this study possible. We also thank N. Simard, S. Laronde, L. Jaimeson, A. Khandani, T. Velanthapillai, and S. Zhao at the UHN/SickKids Flow Cytometry Facility; P. Lo and R. Lopez from the UHN Animal Resource Centre; the Genome Sequencing Facility of the Hartwell Centre for Bioinformatics and Biotechnology, and the Biorepository of St. Jude Children's Research Hospital. We thank J. Ho, E. Lechman, A. Mitchell, L. Jin, M. Doedens, J. Loo-Yong-Kee, K.L. Woo, M.C. Shoong, and J. Roth for their technical assistance; J.A. Kennedy for clinical annotation; and K. Kaufmann for data analysis script. This work was supported by funds to J.E. Dick from the Princess Margaret Cancer Centre Foundation, Ontario Institute for Cancer Research, with funding from the Province of Ontario, Canadian Institutes for Health Research, Canadian Cancer Society Research Institute, Terry Fox Research Institute Program Project, Genome Canada through the Ontario Genomics Institute, and a Canada Research Chair; and to C.G. Mullighan from the the American Lebanese Syrian Associated Charities of St. Jude Children's Research Hospital, a St. Baldrick's Foundation Robert J. Arceci Innovation Award, the Henry Schueler 41&9 Foundation, the NCI grants P30 CA021765 (St. Jude Cancer Center Support Grant) and Outstanding Investigator Award R35 CA197695. Q. Morris was supported by funds from the Natural Sciences and Engineering Research Council, a Compute the Cure award from the NVIDIA charitable foundation, and the Vector Institute. E. Waanders is funded by the Dutch Cancer Society (KUN2012-5366).

1.
Forman
SJ
,
Rowe
JM
. 
The myth of the second remission of acute leukemia in the adult
.
Blood
2013
;
121
:
1077
82
.
2.
Liew
E
,
Atenafu
EG
,
Schimmer
AD
,
Yee
KW
,
Schuh
AC
,
Minden
MD
, et al
Outcomes of adult patients with relapsed acute lymphoblastic leukemia following frontline treatment with a pediatric regimen
.
Leuk Res
2012
;
36
:
1517
20
.
3.
Hunger
SP
,
Mullighan
CG
. 
Acute lymphoblastic leukemia in children
.
N Engl J Med
2015
;
373
:
1541
52
.
4.
Mullighan
CG
,
Phillips
LA
,
Su
X
,
Ma
J
,
Miller
CB
,
Shurtleff
SA
, et al
Genomic analysis of the clonal origins of relapsed acute lymphoblastic leukemia
.
Science
2008
;
322
:
1377
80
.
5.
Ma
X
,
Edmonson
M
,
Yergeau
D
,
Muzny
DM
,
Hampton
OA
,
Rusch
M
, et al
Rise and fall of subclones from diagnosis to relapse in pediatric B-acute lymphoblastic leukaemia
.
Nat Commun
2015
;
6
:
6604
.
6.
Oshima
K
,
Khiabanian
H
,
da Silva-Almeida
AC
,
Tzoneva
G
,
Abate
F
,
Ambesi-Impiombato
A
, et al
Mutational landscape, clonal evolution patterns, and role of RAS mutations in relapsed acute lymphoblastic leukemia
.
Proc Natl Acad Sci U S A
2016
;
113
:
11306
11
.
7.
Waanders
E
,
Gu
Z
,
Dobsom
SM
,
Antic
Z
,
Crawford
JC
,
Ma
X
, et al
Mutational landscape and patterns of clonal evolution in relapsed pediatric acute lymphoblastic leukemia
.
Blood Cancer Discov
2020
;
1
:
1
16
.
8.
Anderson
K
,
Lutz
C
,
van Delft
FW
,
Bateman
CM
,
Guo
Y
,
Colman
SM
, et al
Genetic variegation of clonal architecture and propagating cells in leukaemia
.
Nature
2011
;
469
:
356
61
.
9.
Notta
F
,
Mullighan
CG
,
Wang
JC
,
Poeppl
A
,
Doulatov
S
,
Phillips
LA
, et al
Evolution of human BCR-ABL1 lymphoblastic leukaemia-initiating cells
.
Nature
2011
;
469
:
362
7
.
10.
Kreso
A
,
Dick
JE
. 
Evolution of the cancer stem cell model
.
Cell Stem Cell
2014
;
14
:
275
91
.
11.
Foo
J
,
Michor
F
. 
Evolution of acquired resistance to anti-cancer therapy
.
J Theor Biol
2014
;
355
:
10
20
.
12.
Kim
C
,
Gao
R
,
Sei
E
,
Brandt
R
,
Hartman
J
,
Hatschek
T
, et al
Chemoresistance evolution in triple-negative breast cancer delineated by single-cell sequencing
.
Cell
2018
;
173
:
879
93
.
13.
Iqbal
Z
,
Aleem
A
,
Iqbal
M
,
Naqvi
MI
,
Gill
A
,
Taj
AS
, et al
Sensitive detection of pre-existing BCR-ABL kinase domain mutations in CD34+ cells of newly diagnosed chronic-phase chronic myeloid leukemia patients is associated with imatinib resistance: implications in the post-imatinib era
.
PLoS One
2013
;
8
:
e55717
.
14.
Tzoneva
G
,
Dieck
CL
,
Oshima
K
,
Ambesi-Impiombato
A
,
Sanchez-Martin
M
,
Madubata
CJ
, et al
Clonal evolution mechanisms in NT5C2 mutant-relapsed acute lymphoblastic leukaemia
.
Nature
2018
;
553
:
511
4
.
15.
Hunter
C
,
Smith
R
,
Cahill
DP
,
Stephens
P
,
Stevens
C
,
Teague
J
, et al
A hypermutation phenotype and somatic MSH6 mutations in recurrent human malignant gliomas after alkylator chemotherapy
.
Cancer Res
2006
;
66
:
3987
91
.
16.
Alexandrov
LB
,
Nik-Zainal
S
,
Wedge
DC
,
Aparicio
SA
,
Behjati
S
,
Biankin
AV
, et al
Signatures of mutational processes in human cancer
.
Nature
2013
;
500
:
415
21
.
17.
Churchman
ML
,
Low
J
,
Qu
C
,
Paietta
EM
,
Kasper
LH
,
Chang
Y
, et al
Efficacy of retinoids in IKZF1-mutated BCR-ABL1 acute lymphoblastic leukemia
.
Cancer Cell
2015
;
28
:
343
56
.
18.
Ebinger
S
,
Ozdemir
EZ
,
Ziegenhain
C
,
Tiedt
S
,
Castro Alves
C
,
Grunert
M
, et al
Characterization of rare, dormant, and therapy-resistant cells in acute lymphoblastic leukemia
.
Cancer Cell
2016
;
30
:
849
62
.
19.
Polak
R
,
de Rooij
B
,
Pieters
R
,
den Boer
ML
. 
B-cell precursor acute lymphoblastic leukemia cells use tunneling nanotubes to orchestrate their microenvironment
.
Blood
2015
;
126
:
2404
14
.
20.
Iwamoto
S
,
Mihara
K
,
Downing
JR
,
Pui
CH
,
Campana
D
. 
Mesenchymal cells regulate the response of acute lymphoblastic leukemia cells to asparaginase
.
J Clin Invest
2007
;
117
:
1049
57
.
21.
Mullighan
CG
,
Zhang
J
,
Kasper
LH
,
Lerach
S
,
Payne-Turner
D
,
Phillips
LA
, et al
CREBBP mutations in relapsed acute lymphoblastic leukaemia
.
Nature
2011
;
471
:
235
9
.
22.
Kamel-Reid
S
,
Letarte
M
,
Sirard
C
,
Doedens
M
,
Grunberger
T
,
Fulop
G
, et al
A model of human acute lymphoblastic leukemia in immune-deficient SCID mice
.
Science
1989
;
246
:
1597
600
.
23.
Lock
RB
,
Liem
N
,
Farnsworth
ML
,
Milross
CG
,
Xue
C
,
Tajbakhsh
M
, et al
The nonobese diabetic/severe combined immunodeficient (NOD/SCID) mouse model of childhood acute lymphoblastic leukemia reveals intrinsic differences in biologic characteristics at diagnosis and relapse
.
Blood
2002
;
99
:
4100
8
.
24.
Ng
SW
,
Mitchell
A
,
Kennedy
JA
,
Chen
WC
,
McLeod
J
,
Ibrahimova
N
, et al
A 17-gene stemness score for rapid determination of risk in acute leukaemia
.
Nature
2016
;
540
:
433
7
.
25.
Liem
NL
,
Papa
RA
,
Milross
CG
,
Schmid
MA
,
Tajbakhsh
M
,
Choi
S
, et al
Characterization of childhood acute lymphoblastic leukemia xenograft models for the preclinical evaluation of new therapies
.
Blood
2004
;
103
:
3905
14
.
26.
Meyer
LH
,
Eckhoff
SM
,
Queudeville
M
,
Kraus
JM
,
Giordan
M
,
Stursberg
J
, et al
Early relapse in ALL is identified by time to leukemia in NOD/SCID mice and is characterized by a gene signature involving survival pathways
.
Cancer Cell
2011
;
19
:
206
17
.
27.
Shlush
LI
,
Mitchell
A
,
Heisler
L
,
Abelson
S
,
Ng
SWK
,
Trotman-Grant
A
, et al
Tracing the origins of relapse in acute myeloid leukaemia to stem cells
.
Nature
2017
;
547
:
104
8
.
28.
Clappier
E
,
Gerby
B
,
Sigaux
F
,
Delord
M
,
Touzri
F
,
Hernandez
L
, et al
Clonal selection in xenografted human T cell acute lymphoblastic leukemia recapitulates gain of malignancy at relapse
.
J Exp Med
2011
;
208
:
653
61
.
29.
Klco
JM
,
Spencer
DH
,
Miller
CA
,
Griffith
M
,
Lamprecht
TL
,
O'Laughlin
M
, et al
Functional heterogeneity of genetically defined subclones in acute myeloid leukemia
.
Cancer Cell
2014
;
25
:
379
92
.
30.
Richter-Pechanska
P
,
Kunz
JB
,
Bornhauser
B
,
von Knebel Doeberitz
C
,
Rausch
T
,
Erarslan-Uysal
B
, et al
PDX models recapitulate the genetic and epigenetic landscape of pediatric T-cell leukemia
.
EMBO Mol Med
2018
;
10:e9443
.
31.
Rehe
K
,
Wilson
K
,
Bomken
S
,
Williamson
D
,
Irving
J
,
den Boer
ML
, et al
Acute B lymphoblastic leukaemia-propagating cells are present at high frequency in diverse lymphoblast populations
.
EMBO Mol Med
2013
;
5
:
38
51
.
32.
Schmitz
M
,
Breithaupt
P
,
Scheidegger
N
,
Cario
G
,
Bonapace
L
,
Meissner
B
, et al
Xenografts of highly resistant leukemia recapitulate the clonal composition of the leukemogenic compartment.
Blood
2011
;
118
:
1854
64
.
33.
Chiaretti
S
,
Li
X
,
Gentleman
R
,
Vitale
A
,
Vignetti
M
,
Mandelli
F
, et al
Gene expression profile of adult T-cell acute lymphocytic leukemia identifies distinct subsets of patients with different response to therapy and survival
.
Blood
2004
;
103
:
2771
8
.
34.
Slamova
L
,
Starkova
J
,
Fronkova
E
,
Zaliova
M
,
Reznickova
L
,
van Delft
FW
, et al
CD2-positive B-cell precursor acute lymphoblastic leukemia with an early switch to the monocytic lineage
.
Leukemia
2014
;
28
:
609
20
.
35.
O'Connell
MP
,
Marchbank
K
,
Webster
MR
,
Valiga
AA
,
Kaur
A
,
Vultur
A
, et al
Hypoxia induces phenotypic plasticity and therapy resistance in melanoma via the tyrosine kinase receptors ROR1 and ROR2
.
Cancer Discov
2013
;
3
:
1378
93
.
36.
Kemper
K
,
de Goeje
PL
,
Peeper
DS
,
van Amerongen
R
. 
Phenotype switching: tumor cell plasticity as a resistance mechanism and target for therapy
.
Cancer Res
2014
;
74
:
5937
41
.
37.
Chaidos
A
,
Barnes
CP
,
Cowan
G
,
May
PC
,
Melo
V
,
Hatjiharissi
E
, et al
Clinical drug resistance linked to interconvertible phenotypic and functional states of tumor-propagating cells in multiple myeloma
.
Blood
2013
;
121
:
318
28
.
38.
Pisco
AO
,
Huang
S
. 
Non-genetic cancer cell plasticity and therapy-induced stemness in tumour relapse: ‘What does not kill me strengthens me’
.
Br J Cancer
2015
;
112
:
1725
32
.
39.
Belviso
S
,
Iuliano
R
,
Amato
R
,
Perrotti
N
,
Menniti
M
. 
The human asparaginase enzyme (ASPG) inhibits growth in leukemic cells
.
PLoS One
2017
;
12
:
e0178174
.
40.
Novershtern
N
,
Subramanian
A
,
Lawton
LN
,
Mak
RH
,
Haining
WN
,
McConkey
ME
, et al
Densely interconnected transcriptional circuits control cell states in human hematopoiesis
.
Cell
2011
;
144
:
296
309
.
41.
Jones
CL
,
Stevens
BM
,
D'Alessandro
A
,
Reisz
JA
,
Culp-Hill
R
,
Nemkov
T
, et al
Inhibition of amino acid metabolism selectively targets human leukemia stem cells
.
Cancer Cell
2018
;
34
:
724
40
.
42.
van Galen
P
,
Mbong
N
,
Kreso
A
,
Schoof
EM
,
Wagenblast
E
,
Ng
SWK
, et al
Integrated stress response activity marks stem cells in normal hematopoiesis and leukemia
.
Cell Rep
2018
;
25
:
1109
17
.
43.
Mullighan
CG
,
Su
X
,
Zhang
J
,
Radtke
I
,
Phillips
LA
,
Miller
CB
, et al
Deletion of IKZF1 and prognosis in acute lymphoblastic leukemia
.
N Engl J Med
2009
;
360
:
470
80
.
44.
Luria
SE
,
Delbruck
M
. 
Mutations of bacteria from virus sensitivity to virus resistance
.
Genetics
1943
;
28
:
491
511
.
45.
Eppert
K
,
Takenaka
K
,
Lechman
ER
,
Waldron
L
,
Nilsson
B
,
van Galen
P
, et al
Stem cell gene expression programs influence clinical outcome in human leukemia
.
Nat Med
2011
;
17
:
1086
93
.
46.
Jin
L
,
Hope
KJ
,
Zhai
Q
,
Smadja-Joffe
F
,
Dick
JE
. 
Targeting of CD44 eradicates human acute myeloid leukemic stem cells
.
Nat Med
2006
;
12
:
1167
74
.
47.
Chen
WC
,
Yuan
JS
,
Xing
Y
,
Mitchell
A
,
Mbong
N
,
Popescu
AC
, et al
An integrated analysis of heterogeneous drug responses in acute myeloid leukemia that enables the discovery of predictive biomarkers
.
Cancer Res
2016
;
76
:
1214
24
.
48.
Mazurier
F
,
Doedens
M
,
Gan
OI
,
Dick
JE
. 
Rapid myeloerythroid repopulation after intrafemoral transplantation of NOD-SCID mice reveals a new class of human stem cells
.
Nat Med
2003
;
9
:
959
63
.
49.
Hu
Y
,
Smyth
GK
. 
ELDA: extreme limiting dilution analysis for comparing depleted and enriched populations in stem cell and other assays
.
J Immunol Methods
2009
;
347
:
70
8
.
50.
Edmonson
MN
,
Zhang
J
,
Yan
C
,
Finney
RP
,
Meerzaman
DM
,
Buetow
KH
. 
Bambino: a variant detector and alignment viewer for next-generation sequencing data in the SAM/BAM format
.
Bioinformatics
2011
;
27
:
865
6
.
51.
Mullighan
CG
. 
Single nucleotide polymorphism microarray analysis of genetic alterations in cancer
.
Methods Mol Biol
2011
;
730
:
235
58
.
52.
Pounds
S
,
Cheng
C
,
Mullighan
C
,
Raimondi
SC
,
Shurtleff
S
,
Downing
JR
. 
Reference alignment of SNP microarray signals for copy number analysis of tumors
.
Bioinformatics
2009
;
25
:
315
21
.
53.
Olshen
AB
,
Venkatraman
ES
,
Lucito
R
,
Wigler
M
. 
Circular binary segmentation for the analysis of array-based DNA copy number data
.
Biostatistics
2004
;
5
:
557
72
.
54.
Venkatraman
ES
,
Olshen
AB
. 
A faster circular binary segmentation algorithm for the analysis of array CGH data
.
Bioinformatics
2007
;
23
:
657
63
.
55.
Lin
M
,
Wei
LJ
,
Sellers
WR
,
Lieberfarb
M
,
Wong
WH
,
Li
C
. 
dChipSNP: significance curve and clustering of SNP-array-based loss-of-heterozygosity data
.
Bioinformatics
2004
;
20
:
1233
40
.
56.
Dobin
A
,
Davis
CA
,
Schlesinger
F
,
Drenkow
J
,
Zaleski
C
,
Jha
S
, et al
STAR: ultrafast universal RNA-seq aligner
.
Bioinformatics
2013
;
29
:
15
21
.
57.
Anders
S
,
Pyl
PT
,
Huber
W
. 
HTSeq–a Python framework to work with high-throughput sequencing data
.
Bioinformatics
2015
;
31
:
166
9
.
58.
Love
MI
,
Huber
W
,
Anders
S
. 
Moderated estimation of fold change and dispersion for RNA-seq data with DESeq2
.
Genome Biol
2014
;
15
:
550
.
59.
Reimand
J
,
Isserlin
R
,
Voisin
V
,
Kucera
M
,
Tannus-Lopes
C
,
Rostamianfar
A
, et al
Pathway enrichment analysis and visualization of omics data using g:Profiler, GSEA, Cytoscape and EnrichmentMap
.
Nat Protoc
2019
;
14
:
482
517
.
60.
Notta
F
,
Doulatov
S
,
Laurenti
E
,
Poeppl
A
,
Jurisica
I
,
Dick
JE
. 
Isolation of single human hematopoietic stem cells capable of long-term multilineage engraftment
.
Science
2011
;
333
:
218
21
.