AALL0434 participant identification and clinical annotation

COG studies AALL0434 (NCT00408005) and AALL1231 (NCT02112916) were approved by the National Cancer Institute Cancer Evaluation and Therapeutic Program, the US Food and Drug Administration, the Pediatric Central Institutional Review Board (IRB) and local IRBs at all participating centers. Written informed consent was obtained from study participants and, when appropriate, their legal guardians, in accordance with the Declaration of Helsinki. All participant data were deidentified and written informed consent was obtained to publish the indirect identifiers in the present manuscript. Further information on research design is available in the Nature Portfolio Reporting Summary linked to this article.

Secondary genomic studies were approved by the Children’s Hospital of Philadelphia (CHOP) IRB. In total, 40 cases from AALL0434 (Supplementary Table 1) and eight healthy thymus and bone marrow controls (Supplementary Tables 5 and 6) were selected for single-cell study. The healthy thymus and bone marrow used for this work were residual tissues after collection for clinical care. Leukemia samples were bone marrow or blood samples collected and banked for COG trials. Within AALL0434, ETP status was centrally assessed in diagnostic bone marrow or peripheral blood samples using 8–9-color multiparameter flow cytometry26. ETP was defined as having lymphoblasts that were CD8−CD1a− (<5% positive), weakly expressed CD5 (either <75% positive or median intensity more than one log less than mature T cells) and expressed one or more myeloid or stem cell markers (>25% positive) including CD13, CD33, CD34, CD117 and HLA-DR (ref. 13). Subjects meeting the ETP immunophenotypic criteria but with stronger expression of CD5 were classified as near-ETP. Subjects with neither ETP nor near-ETP were defined as non-ETP. MRD was assessed using 8–9-color flow cytometry and was performed using established methods at a COG flow cytometry reference laboratory (University of Washington or Johns Hopkins University).

Processing of T-ALL diagnosis samples

Peripheral blood or bone marrow aspirate samples were thawed at 37 °C, treated with 1:10 (v/v) 1 mg ml−1 DNase I (Sigma-Aldrich, D4513) for 90 s at 37 °C, resuspended in 10 ml of Iscove’s modified Dulbecco’s medium (IMDM) + 2% FBS and centrifuged (160g for 5 min). Samples were retreated with DNase I and resuspended in fluorescence-activated cell sorting (FACS) buffer (Ca2+-free and Mg2+-free PBS + 1% BSA). Cell number and viability were recorded using a Countess II cell counter (Invitrogen). More than 1 million live cells were aliquoted for tail-vein injection into NSG mice, with the remaining stained with DAPI (Invitrogen, D1306) and subjected to FACS sorting (FACSAria Fusion, BD).

scRNA-seq and CITE-seq library preparation

FACS-sorted DAPI-negative live cells were centrifuged and resuspended in cell staining buffer (BioLegend, 420201) at 45 μl per million cells. Cells were blocked with Human TruStain FcX (BioLegend, 422301) at 5 μl per million cells (4 °C, 15 min). After blocking, cells were stained with a TotalSeq-A antibody cocktail (30 min, 4 °C). Cells were washed three times using cell staining buffer (BioLegend, 420201) and resuspended in PBS + 0.04% BSA. Cells were counted using a Countess II cell counter. A total of 17,000 cells per sample were then loaded onto 10x Genomics Chromium controller and processed with the Chromium NEXT GEM single-cell 3′ reagent kit (version 3.1). GEX libraries were constructed using the 10x Genomics library preparation kit following the instructions. Antibody-derived tag (ADT) libraries were constructed using the KAPA HiFi HotStart ReadyMix kit (Kapa Biosystems, KK2601). The following program was used for ADT library PCR: 98 °C for 2 min, 14–15 cycles of 98 °C for 20 s, 60 °C for 30 s and 72 °C for 20 s, followed by 72 °C for 5 min and a hold at 4 °C. Library quality was checked using the Agilent high-sensitivity DNA kit (Agilent, 5067-4626) and Bioanalyzer 2100. Libraries were quantified using the dsDNA high-sensitivity assay kit (Invitrogen, Q33231) on a Qubit fluorometer and quantified using the qPCR-based KAPA quantification kit (Kapa Biosystems, KK4844). Libraries were sequenced on an Illumina NovaSeq 6000 with 28:8:0:87 paired-end format.

scATAC-seq library preparation

DAPI-negative live cells were centrifuged at 300g (5 min at 4 °C), mixed in 45 μl of lysis buffer and incubated (3 min on ice). Next, 50 μl of prechilled wash buffer was added without mixing and centrifuged immediately at 300g (5 min at 4 °C). Then, 95 μl of supernatant was discarded, 45 μl of diluted nuclei buffer (10x Genomics) was added and the sample was centrifuged (300g; 5 min at 4 °C). The nuclear pellet was then resuspended in 7 μl of prechilled diluted nuclei buffer and the nuclear concentration was determined using a Countess II cell counter. A total of 7,000–20,000 nuclei were used for the transposition reaction in bulk, loaded onto the 10x Genomics Chromium controller and processed with the Chromium NEXT GEM scATAC reagent kit (version 1.1). Library quality was checked using the Agilent high-sensitivity DNA kit and Bioanalyzer 2100. Libraries were quantified using the dsDNA high-sensitivity assay kit on a Qubit fluorometer and quantified using the qPCR-based KAPA quantification kit. Libraries were sequenced on an Illumina NovaSeq 6000 with 49:8:16:49 paired-end format.

Expansion and profiling of T-ALL blasts in PDX

NSG mice (RRID:IMSR_JSX:005557) were used for all experiments. For the development of PDX models, we injected ~106 blasts from viably frozen participant samples (bone marrow or blood) per mouse to develop primagrafts (Supplemental Table 22). PDX-expanded blasts were isolated from the spleen or bone marrow. Frozen samples were thawed (37 °C), resuspended in IMDM + 2% FBS and treated with DNase I twice. Cells were washed twice with RPMI medium, resuspended in flow buffer, stained with DAPI and anti-human CD45 antibody (BD Pharmingen, 555485) and subjected to FACS sorting (FACSAria Fusion, BD). DAPI-negative hCD45+ sorted cells were stained with 10x Genomics 3′ CellPlex multiplexing solution, washed three times and immediately processed using the 10x Genomics Chromium controller and the Chromium NEXT GEM single-cell 3′ reagent kit (version 3.1). The 3′ GEX libraries were constructed using the 10x Genomics library preparation kit. CellPlex libraries were constructed using the 10x Genomics 3′ CellPlex kit. Library quality was checked using the Agilent high-sensitivity DNA kit and Bioanalyzer 2100. Libraries were quantified using the dsDNA high-sensitivity assay kit on a Qubit fluorometer and quantified using the qPCR-based KAPA quantification kit. Libraries were sequenced on an Illumina NovaSeq 6000 with 28:8:0:87 paired-end format.

CD34+ progenitor isolation from infant or pediatric thymi

Pediatric thymi were obtained and used according to and with the approval of the CHOP IRB. Thymus tissue was mechanically disrupted and treated with liberase (0.2 mg ml−1, 30 min at 37 °C; Roche) with intermittent shaking, as previously described19. Thymocytes were resuspended into flow buffer, sorted into DAPI-negative lineage-negative CD34+CD1A− fractions and subjected to scRNA-seq and scATAC-seq.

Projection onto healthy reference trajectory

Participant-derived cells were projected onto the healthy reference trajectory using the MapQuery function in Seurat 4.0.5. For scRNA-seq data, participant and healthy control data were coembedded into a low-dimensional space using the default anchor-based canonical correlation analysis (CCA) method in Seurat 4.0.5 (30 dimensions, 2,000 anchor features) and cell type label transfer was performed on a sample-by-sample basis using the TransferData function. For scATAC-seq data, peaks from participant and healthy reference data were merged using the mergePeaks module from scATAC-pro53 and peak × cell matrices with merged peaks were reconstructed for each participant with the scATAC-pro reConstMtx module. This allowed for participant and healthy control data to be coembedded into a low-dimensional space analogous to the scRNA-seq data.

AALL0434 ETP-ALL stratification using BMP-119 signature

BMP-like and T-specified DEGs were stringently filtered using cutoffs of false discovery rate (FDR) < 0.001 and average log2 fold change (log2FC) > 0.5, leaving 66 BMP-like DEGs and 53 T-specified DEGs. The z-score-based signature scoring was performed on 110 bulk-sequenced participants with ETP-ALL with BMP-like DEGs as positive features and T-specified DEGs as negative features. For each participant, the mean T-specified feature z-score was subtracted from the mean BMP-like feature z-score, with a score of >0 being interpreted as more BMP-like than T-specified. This cutoff was selected to compare participants on the basis of a relative enrichment of either phenotype. Participants were binarized by BMP-like signature score (BMP-like > T-specified versus T-specified > BMP-like) and OS and EFS were compared using the Cox proportional hazard model with day 29 MRD and CNS status taken as covariates using the survfit function from survival 3.2-13 (‘survfit(Surv(time.survival, status.survival) ~ high.BMP + D29.MRD + D29.CNS.status’).

Integration of single-cell signatures with mutation calls

Bulk RNA-seq data for n = 110 ETP samples with corresponding WES and WGS mutation calls were scored using 66 BMP-like DEGs and 53 T-specified DEGs using AUCell 1.12.0. For 1,490 mutant genes in 110 ETP samples, the number of samples carrying mutations was quantified and the mean BMP-like area under the curve (AUC) and T-specified AUC were calculated. Mutant genes observed in ≥5 samples with mean VAF > 0.05 were plotted for visualization. Classification of genes was derived from a previous bulk genomics study on ETP-ALL. For fusion drivers, the mean BMP-like AUC, the mean T-specified AUC, the percentage of participants with positive EOI MRD, the percentage of participants that died during the trial and the number of unique fusion partners were calculated.

Identification of a consensus BMP-like gene signature

BMP-like DEGs from participants with ETP-ALL (n = 56 BMP-like versus T-specified) and participants with non-ETP-ALL (n = 445 BMP-like versus postcommitment) were overlapped and the average log2FC was calculated. A total of 17 genes with average log2FC > 0.9 were retained as a consensus BMP-17 signature. We performed AUC-based signature scoring using AUCell 1.12.0 (with the top 25% of expressed genes) on bulk RNA-seq diagnostic T-ALL samples from two independent COG trials using BMP-17 DEGs. We then binarized participants on the basis of AUCell signature score and used the Cox proportional hazard model with EOI MRD and CNS status taken as covariates using the survfit function from survival 3.2-13 (‘survfit(Surv(time.OS, status.OS) ~ BMP-17 + high.BMP-17 + D29.MRD + CNS.status). In each case, the top half of participants was compared to the bottom half of participants.

In silico drug screening against BMP-like blasts

Drug–target data from two independent drug target databases (TTD54 and DrugIDB55) and a third database (OpenTargets56) that focuses on next-generation targets were overlapped with BMP-like DEGs (log2FC > 0.2; adjusted P < 0.01). Targetable gene products were given a score of 1 for each database in which a resulting hit was obtained. To search for drugs that could specifically modify the BMP-like state, we inputted top BMP-like DEGs and TFs (n = 56) and top T-specified DEGs and TFs (n = 62) into the LINCS1000 (ref. 57) database under default parameters. Perturbation results were filtered in R to filter compound-mediated perturbations for compounds with defined targets, statistical significance (log10FDR > 1), effect size (normalized connectivity score > 0.8), specificity to BMP-like state (raw connectivity score > 0) and activity in two or more leukemia cell lines. Non-compound perturbations were filtered for statistical significance (log10FDR > 1) and effect size (normalized connectivity score > 0.8) and further separated into gene overexpression and gene knockdown (including short hairpin RNA knockdown, clustered regularly interspaced short palindromic repeats knockout and ligand-based perturbation) classes. BMP-like DEGs targeted by top compound perturbations and/or genes with overexpression or knockdown were given a score of 1. Lastly, we identified BMP-like DEGs that showed increased dependency in leukemia cell lines (n = 59) compared to non-leukemia and non-lymphoma cell lines (n = 1,052) in the cancer dependency map58 (DepMap) portal. Genes with negative dependency scores in leukemia cell lines (mean dependency score < −0.1), dependency FC > 2 and >25% expression in BMP-like blasts were assigned a score of 1. Next, BMP-like DEGs with log2FC > 1 were assigned a score of 1 and genes with log2FC between 0.5 and 1 were assigned a score of 0.5. Beyond high expression change, we prioritized BMP-like DEGs with high percentage expression in BMP-like blasts; genes with >80% expression were assigned a value of 1, whereas genes with 50–80% expression were assigned a value of 0.5. Finally, we ranked genes with high statistical significance (adjusted P < 1 × 10−100 was given a score of 1; adjusted P < 1 × 10−50 was given a score of 0.5). The sum of DE evidence (ranging from 0 to 3) and database evidence (ranging from 0 to 5) was taken to rank BMP-like DEGs for follow-up experimental studies.

In vitro drug screening with leukemia active drug panel

Human leukemia blasts were collected from mouse spleen and enriched using an immunomagnetic isolation kit (StemCell Technologies, 19849) and screened with a panel of 40 leukemia active drugs (Supplementary Table 18) using a previously described imaging-based assay with a stromal cell coculture system59.

Nomination of BMP-like specific drugs from drug screening

PDX-expanded blasts from five participants were screened in a stromal coculture system and dose–response curves were generated for each compound with the primary readout being cell viability (as a percentage of control). We defined ETP active drugs as compounds with IC50 < 1,000 nM and categorized each compound as not active (active in 0/4 participants with ETP), partially active (active in 1–3 participants with ETP) or active (active in 4/4 participants with ETP). We then compared IC50 values for ETP active compounds and used three comparisons to nominate drugs that were differentially active in high-BMP participants: BMP-high and MRD-positive (n = 3) versus BMP-low and MRD-negative (n = 2: one ETP and one non-ETP); BMP-high and MRD-positive (n = 2) versus BMP-low (n = 1); BMP-high and MRD-positive (n = 3) versus BMP-low (n = 1). Drugs with differential activity in BMP-high participants in all three comparisons were nominated as BMP-specific drugs. The sensitivity of these drugs was confirmed using PDX-expanded blasts from five additional participants, for a total of five BMP-high models and five BMP-low models.

scRNA-seq and CITE-seq data processing

Demultiplexing and alignment of RNA and ADT sequences were performed with Cell Ranger 3.1.0. Low-quality cells and red blood cells were then filtered by retaining only cells with between 300 and 2,500 genes in the scRNA-seq data, greater than 1,500 RNA counts, less than 10% mitochondrial RNA and fewer than three unique molecular identifiers (UMIs) mapping to hemoglobin B. To remove cell doublets in scRNA-seq data, DoubletFinder 2.0.3 was run with 5% of the expected rate of doublets. Participant cell × gene and cell × ADT count matrices were individually saved and subsequently concatenated using Seurat 4.0.5 for downstream analyses. For some analyses as specified, cell × gene and cell × ADT matrices for each participant were subset for G1 cells (representing the phenotype most resistant to conventional therapy) and were downsampled to match the lowest value in the cohort. After log-normalization, the FindVariableFeatures function in Seurat 4.0.5 was used to identify the top 5,000 features with greater than expected variance. Variable features with expression in >1% of cells were kept as the input to principal component analysis (PCA), with subsequent visualization being performed using uniform manifold approximation and projection (UMAP) of the top 50 principal components (PCs), 30 neighbors and two PCs. For visualization, we used the IntegrateLayers function in Seurat 5.0.3 with the RPCAIntegration method and default parameters.

scATAC-seq data processing

Demultiplexing of scATAC-seq reads was performed with Cell Ranger-ATAC 1.1.0 (alignment to hg38) and peak calling was performed with BWA and MACS2 using the scATAC-pro pipeline53 with default parameters. Low-quality cells were filtered for those cells with <3,000 (low quality) and >50,000 unique fragments (doublets), <40% reads in peaks (fraction of reads in peaks < 0.4) and >20% reads mapping to mitochondria. To construct a common peak set, the top 100,000 peaks (defined by MACS2 MapQ score) were selected for downstream merging, alongside 1,500 randomly selected cell barcodes from each participant. We defined two sample sets for merging peaks: one with 25 participants with ETP-ALL and one with 40 participants with T-ALL. For each sample set, peaks were merged with the scATAC-pro mergePeaks module and peak × cell matrices with merged peaks and downsampled cell barcodes were reconstructed with the scATAC-pro reConstMtx module.

Bulk RNA-seq analysis and visualization

Sequencing read adaptors were removed using Trim Galore 0.4.4 (https://www.bioinformatics.babraham.ac.uk/projects/trim_galore/) with parameters ‘-q 20 --phred 33 --paired’. Reads were aligned to the human genome GRCh38 using STAR. The resulting gene counts for each sample were estimated by RSEM, and combined as gene count matrix. RSEM expected counts were processed and filtered. First, samples were required to exhibit expression of over one count per million in ≥5 samples. Second, batch correction was performed using the sva R package 3.46 function ComBat_seq. Batches were defined on the basis of library type as stranded or unstranded and on the basis of cohort as TARGET or X01 sequenced samples. Third, the DESeq2 R package 1.38.3 vst function was used for data normalization. Limma 3.46.0 was used for DE analysis. For visualization, the raw counts were transformed into transcripts per million and visualized using UMAP with the top ten PCs and k = 30 neighbors.

Construction of healthy reference trajectory

Construction of the healthy reference trajectory began with sample-by-sample cell annotation followed by consensus clustering and annotation. Annotations from previously published bone marrow samples were kept60. Cell × gene matrices from each thymus sample were log-normalized and subjected to dimensionality reduction. Cells were clustered at multiple resolutions (k = 1, 2 and 3) and clusters were given preliminary labels on the basis of marker gene expression. Cell × gene matrices from all thymus donors (n = 3) were then concatenated, log-normalized and subjected to dimensionality reduction. Cells were reclustered at high resolutions (k = 3) and clusters were relabeled on the basis of marker gene expression and prior labels.

scRNA-seq

Cell × gene matrices from healthy thymus donors were then concatenated with cell × gene matrices from healthy bone marrow donors (n = 5), log-normalized and subjected to dimensionality reduction using the top 25 PCs. The FindVariableFeatures function in Seurat was used to identify top 2,000 variable genes. These 2,000 genes were then filtered in two iterations on the basis of the Gini coefficient61. Briefly, a shared nearest neighbor graph was constructed using 50 and 20 PCs; cells were clustered at k = 0.1 resolution and the Gini coefficient was calculated for each variable gene. Genes with a low Gini coefficient (bottom 10% percentile) and cluster level expression < 10% were removed in each iteration. The 134 cell-cycle-related genes previously described19 were removed. The remaining 931 variable features were used as input to PCA and UMAP dimension reductions (25 PCs). Trajectory analysis was performed using Slingshot 1.8.0 with HSPC as the start cluster and effector T, mature B and monocyte as the end clusters for T, B and myeloid trajectories, respectively. Principal curves were selected for T cell and myeloid cell trajectories and values were scaled to a maximum of 1 in each curve. Pseudotime values of shared cell states that occurred in both myeloid and T cell development (multipotent progenitors: HSPC and LMPP) were then averaged. Pseudotime values were scaled for myeloid development (0 to −1) and T cell development (0 to 1). Statistical comparisons in the overall arrest state were made using a two-sample Kolmogorov–Smirnov test, as previously described62.

scATAC-seq

Gene–activity matrices for scATAC-seq were constructed by summing counts within the gene body and 2 kb upstream, as previously described60. Integration of scATAC-seq samples was performed using gene–activity matrices and Seurat 4.0.5 using the default anchor-based CCA method with 30 dimensions, 2,000 anchor features and k.filter = 100. To learn labels for scATAC-seq data from scRNA-seq data, transfer anchors were computed using CCA with scRNA-seq as the reference and cell type label transfer was performed on a sample-by-sample basis using the TransferData function. Cell × peak matrices from all thymus donors were then concatenated. Cells were reclustered at high resolution (k = 3) and clusters were reannotated according to consensus labels. Dimensionality reduction was performed using UMAP of the top ten PCs of the concatenated scATAC-seq data and trajectory analysis was performed as described above.

Level 1 annotation of CITE-seq data

To distinguish malignant blasts from non-malignant cells, we first used a cluster-based statistic, Shannon entropy, to identify clusters of cells at multiple clustering resolutions (k = 1, 2 and 3) to identify four cell populations that had contribution from every participant. The Shannon entropy statistic was calculated using the formula \(-\sum }\!\left(}\right)\times \log }\!\left(}\right)\), where p(x) is the frequency of cells arising from any one participant in any one cluster, ranging from 0 to 1. Second, we concatenated and clustered participant-derived single-cell data with healthy bone marrow and thymus controls. Third, we calculated a similarity score to healthy controls across all participant-derived cells. Participant data and healthy control data were coembedded into a low-dimensional space using the default anchor-based CCA method in Seurat 4.0.5 (30 dimensions and 2,000 anchor features) and a k = 30 mutual nearest neighbor score was assigned for each cell using the TransferData function. Copy number profiles were analyzed using InferCNV 1.6.0 on a randomly downsampled (1:10) subset of participant data. We then compared blast percentages calculated in scRNA-seq to pathology reports of blast percentage obtained from diagnostic aspirate (mean absolute deviation = 8.9%; non-significant difference according to paired two-tailed t-test).

Level 1 annotation of single-nucleus (sn)ATAC-seq data

Firstly, annotated scRNA-seq data were used as a reference to annotate participant scATAC-seq data on a paired, sample-by-sample basis. For each participant, gene–activity matrices for scATAC-seq were constructed by summing counts within the gene body and 2 kb upstream, as previously described63. Integration of scATAC-seq samples with scRNA-seq data was performed using gene–activity matrices and Seurat 4.2.0 using the default anchor-based CCA method with 30 dimensions, 2,000 anchor features and k.filter = 100 using the TransferData function on a sample-by-sample basis. Then, participant data and healthy control data were coembedded into a low-dimensional space using the default anchor-based CCA method in Seurat 4.0.5 (30 dimensions, 2,000 anchor features and k.filter = 100) and a k = 30 mutual nearest neighbor score was assigned using the TransferData function to assess their similarity. Lastly, blast percentages calculated in scATAC-seq were compared to blast percentages calculated in scRNA-seq, showing high concordance (median deviation = 1.2%; non-significant difference according to paired two-tailed t-test).

Differential activity analyses

For TF motif enrichment analysis, cell × deviation score matrices were generated using the addGCBias, matchMotifs, getBackgroundPeaks and computeDeviations functions in chromVAR 1.12.0 with hg38 as the reference genome. Differential activity analysis was performed using the Wilcoxon rank-sum test with Benjamini–Hochberg multiple-testing correction with downsampling. For each motif in any particular comparison, we calculated the Δ median chromVAR deviation score, Δ mean chromVAR deviation score, adjusted P value, percentage expression of corresponding TF in paired scRNA-seq data and the ratio of median and mean chromVAR deviation score. DA TF motifs were defined by Δ median chromVAR deviation score > 0.0025, FDR < 0.001, >20% cell expression of corresponding TF and a ratio of median and mean chromVAR deviation score between 0.7 and 1.3, unless otherwise specified.

Subtype-specific transcriptional regulatory analysis

As we previously described60, for each cell in the scRNA-seq dataset, an scRNA-seq and scATAC-seq ‘metacell’ was defined by pooling counts for each gene or peak from the 30 nearest neighbors in the PC space by cosine distance. Metacell counts were log-normalized and scaled. For a gene of interest, we ran a linear regression model using metacell gene expression as the dependent variable and putative enhancer peaks within 500 kb of the transcription start site as regressors. Bonferroni-adjusted P values < 0.01 with a positive coefficient were considered significant. Top induced targets of TCF7 and LEF1 were defined by high-confidence EP regression (regression coefficient > 0.3) and log2FC > 0.5 for z-scoring on bulk RNA-seq data.

Promoter–enhancer coaccessibility networks (CCANs)

Cicero64, which identifies coaccessible pairs of DNA elements, was implemented in Signac65 through the make_cicero_cds function followed by the run_cicero function with the following parameters: sample_num = 100, window = 500,000. These links were aggregated into cis-coaccessible networks using the generate_ccans function with default parameters. The BMP-like and T-specified CCANs were isolated by identifying links to regions within 2,000 bp of the transcription start sites for the 66 and 53 DEGs for each state, respectively. Any peaks that overlapped regions within that coaccessibility group were then subset as potential regulators. This yielded 660 peaks in the BMP-like CCAN and 1,011 peaks in the T-specified CCAN. Those peaks were then used as input to HOMER with parameters ‘-size 200 -mask’ to identify motifs enriched in coaccessible regions. Motifs with a q value < 0.05 were considered significant.

Transcriptional regulatory analysis of developmental states

An integrated enhancer-driven transcriptional regulatory analysis was conducted using SCENIC+ 1.0a1 (ref. 44) following the standard vignettes with minor modifications. BMP-like and T-specified states from the scRNA-seq data and scATAC-seq data were extracted and 35 topics were empirically selected. The SCENIC+ pipeline was then run in non-multiome mode, using five cells per metacell. The search space was defined as 0–500 kb. Regulons were filtered with the following parameters: rho_threshold = 0.03, min_regions_per_gene = 0 and min_target_genes = 10. All other parameters were maintained as the defaults. Region-based and gene-based specificity scores for the BMP-like and T-specified states were calculated using the regulon_specificity_scores function.

AUCell pathway analysis and GSEA

Pathway analysis was conducted using two methods. First, pathway enrichment scores for gene signatures were defined from our single-cell analysis, including the gene sets defined through our transcriptional regulatory analyses or downloaded from the Molecular Signatures Database66, and were determined using AUCell 1.12.0 with the top 5% of genes. Additional gene sets for NOTCH activation were previously published67,68. GSEA was conducted to compare DE pathways between blast populations. A full gene list was constructed using the FindMarkers function in Seurat 5.0.3 with the following parameters: min.pct = 0.001, logfc.threshold = 0, only.pos = FALSE, max.cells.per.ident = 1,500. This gene list was sorted by log2FC to use as input to preranked GSEA using the fgsea package69.

Cell-cycle analyses in single-cell-sequenced participants with ETP-ALL

Cell-cycle signature scoring and phase classification was performed on ten high-MRD and ten MRD-negative participants with ETP using the CellCycleScoring function in Seurat 4.0.5 with default parameters. A total of 43 S-phase and 54 G2M-phase signature genes70 were provided as input. Cells were then randomly downsampled so that each participant would be represented by an equal cell number (3,350 per participant and 33,500 per group).

BMP-like DE analyses

BMP-like and T-specified DEGs were computed using the FindMarkers function in Seurat 4.0.5 with the following parameters: assay = RNA, logfc.threshold = 0, ident.1 = T-specified-R (T-specified blasts from ten MRD-negative participants), ident.2 = BMP-like-NR (BMP-like blasts from 15 MRD-positive participants) and max.cells.per.ident = 1,500. The input matrix to DE analysis was a matrix of G1-phase ETP-ALL blasts with an equal number of cells per participant (1,711 per participant and 42,775 cells total). To identify DE TFs and DE surface markers, the same process was repeated with using genes encoding human TFs71,72 (feature = TFs) and a change of assay to normalized ADT count matrix (assay = ADT).

Intersection of DE TF and DA motifs

DE TFs from scRNA-seq data, defined by average log2FC > 0.15 and FDR < 0.001, were intersected with DA TF motifs from scATAC-seq data. DA TF motifs were defined by Δ median chromVAR deviation score > 0.0025, FDR < 0.001, >20% cell expression of corresponding TF and a ratio of median and mean chromVAR deviation score between 0.7 and 1.3.

Identification of NOTCH1 mutations in scRNA-seq

Samples were demultiplexed into FASTQ files using bcl2fastq. FASTQ files were then processed using IronThrone 2.1 with the default parameters and inputs for 10x version 3.1 scRNA-seq data. Specifically, for each variant, IronThrone was run in circularization mode (--run = circ) with UMI length 12 (--umilen 12) and cell barcodes from each sample’s Cell Ranger output (--whitelist sample.specific.barcodes.tsv), following the configuration set within IronThrone 2.1 documentation (https://github.com/dan-landau/IronThrone-GoT).

Identification of BMP-like blasts in participants without ETP

Precommitted blasts in ten participants without ETP (six EOI-MRD-negative participants, 7,152 precommitted blasts; four EOI-MRD-positive participants, 11,047 precommitted blasts) were subsetted (total of 52,971 blasts and 15,830 blasts, respectively) and mean proportions for corresponding cell fractions (BMP-like, MEP-like and pro-T cell-like) were quantified in each participant. The mean proportion of each cell type for participants of each group was plotted, with the proportion of BMP-like blasts in MRD-negative versus MRD-positive participants being compared using the prop.test function.

Single-cell signature-based stratification of non-ETP cases

Precommitment and postcommitment DEGs were computed as described above. The input matrix to DE analysis was a matrix of G1-phase non-ETP-ALL blasts with a maximum of n = 1,500 cells per participant (34,384 cells in total). DEGs located on the X and Y chromosomes were filtered out to retain the core biology of both cell fractions. A z-score-based signature scoring was performed on 1,051 bulk-sequenced diagnostic participants with non-ETP-ALL with BMP-like DEGs as positive features and T-specified DEGs as negative features. Human TFs were previously curated71,72. Survival was analyzed using Cox proportional hazards as described above.

Identification of a BMP-like surface marker signature

DE ADTs from ETP-ALL BMP-like blasts and non-ETP BMP-like blasts were overlapped and the average log2FC was calculated. Nine genes with |log2FC| > 0.5 (five with positive expression and four with negative expression) were retained as a consensus BMP-surface-9 signature. Gene signature scoring and survival analysis were conducted as described above.

LASSO (least absolute shrinkage and selection operator) optimization of prognostic gene signatures

Refined gene sets were found using LASSO penalized regression. The Cox proportional hazards model was used with gene z-scores as features (glmnet and survival R packages). The model was adjusted for participant age, white blood cell (WBC) count at diagnosis, CNS status and treatment protocol, by including these as covariates on which no penalty was applied. The penalty was only applied to the gene features but the range of the predicted coefficients was bounded such that genes enriched in the BMP signature were given positive coefficients (increased hazard) and genes enriched in the T-specified signature were given negative coefficients (decreased hazard). The models were fit to the entire AALL0434 RNA-seq dataset and ETP status was used to stratify the survival allowing for different baseline hazards, followed by 100-fold cross-validation to determine the penalty parameter with the lowest root-mean-squared error.

Identification of BMP-like signature expression patterns

To contextualize BMP-17 and BMP-surface-9 marker genes within normal hematopoiesis, we performed AUC-based signature scoring of healthy donor scRNA-seq reference maps using AUCell with the top 10% of expressed genes considered for computational efficiency (aucMaxRank = 0.1). BMP-surface-9 marker genes were divided into positive DEGs (n = 4) and negative DEGs (n = 5) and AUCs were calculated for each gene set. Overall enrichment (that is, the aggregate AUC) of the BMP-surface-9 signature was calculated by taking the difference in AUC between positive and negative surface markers.

Integration of AALL0434 and Lee et al.’s bulk RNA-seq data

In vitro drug sensitivity data were integrated from the current study (ten PDX models) with previously published drug screening data47. To generate a consensus BMP-like signature across data, we first converted the AALL0434 bulk transcriptomic data to fragments per kilobase of transcript per million mapped reads format using the convertCounts function in the DGEobj.utils package. A BMP-like gene signature was scored using the 119 BMP-like and T-specified DEGs with a robust z-score. The negative features (T-specified DEGs) were subtracted from the positive features (BMP-like DEGs).

Statistics and reproducibility

No statistical method was used to predetermine sample size. All data meeting the quality control threshold were included. The investigators were not blinded to allocation during genomic profiling and assessment of participant data. No data points were excluded from analyses related to single-cell and bulk genomics. No animal models were excluded from PDX-related analyses. Randomization and blinding were used for all in vitro and in vivo experiments. Statistical comparisons were made using a two-sided Wilcoxon rank-sum test unless otherwise specified in the figure legend. As the Wilcoxon test is non-parametric, we did not formally test for normality of the data. The chi-squared test was used to compare cell type proportions. The Cox proportional hazards model was used for the assessment of survival outcomes.

Reporting summary

Further information on research design is available in the Nature Portfolio Reporting Summary linked to this article.