Our study includes a cohort of six pediatric patients diagnosed with high-risk B-ALL who underwent early intensification treatment using the CAML protocol (Supplementary Table S1). High-risk B-ALL was diagnosed based on current Chinese pediatric ALL risk stratification guidelines. Peripheral blood samples were collected post-treatment, with a focus on the second round of treatment where three patients achieved remission and three experienced relapses. Additionally, peripheral blood samples from three age-matched healthy controls were obtained for comparative analysis.

Peripheral blood mononuclear cells (PBMCs) were isolated from the collected samples using density gradient centrifugation at 600×g for 25 minutes at 20°C. Cell number and viability were assessed by fluorescence using Reward C100 cell counter. Informed consent was obtained from all participants or their guardians, and the study was approved by the Institutional Review Board.

PBMCs were subjected to scRNA-seq to capture the transcriptomic landscape of individual cells. Single-cell RNA sequencing libraries were prepared using the Chromium Single Cell 3’V2 Reagent Kit, and generated using the 10x Genomic platform. All libraries were sequenced on Illumina NovaSeq 6000 platform. Each sample generated approximately 110 – 140 million raw bases, with Q30 scores above 91%. Quality control, data normalization, and dimensionality reduction were performed using established computational methods.

Upon receiving the FASTQ files from the 10x Genomics single-cell RNA sequencing results, we utilized Cell Ranger software (v7.1.0, 10x Genomics, USA) to align barcodes and UMIs to the reference genome (GRCh38), followed by filtering and quantifying gene expression at the single-cell level, resulting in a cell-gene expression matrix. Putative doublets were identified and removed using Scrublet (Scanpy pipeline) and DoubletFinder (Seurat workflow), based on the expected doublet rate from 10x Genomics specifications. This matrix was then processed using Seurat (v4.3.0, Satija Lab, USA) and Scanpy (v1.9.3, Theis Lab, Germany) (6). Cells were retained if they exhibited between 200 and 5,000 detected genes, 500 to 25,000 unique molecular identifiers (UMIs), and less than 10% mitochondrial gene expression. To further ensure data quality, doublets were identified using DoubletFinder (v2.0.3), with an expected doublet rate of 8% based on cell loading density and a pK value of 0.09. Cells with a doublet prediction score greater than 0.25 were excluded prior to downstream normalization. Cells with fewer than 200 detected genes or with mitochondrial gene read counts exceeding 10% were excluded. We integrated the data across different samples using the IntegrateLayers function through Canonical Correlation Analysis (CCA). Dimensionality reduction was performed via Principal Component Analysis (PCA), and the resulting data were clustered using UMAP. The clusters were annotated using ScType (v1.1.0, Ianevski Lab, Finland) (7), and the expression of marker genes for various immune cell types was examined. Following annotation, the R software (v4.3.1, R Foundation for Statistical Computing, Austria) was employed to calculate the proportion of each cell type across sample groups.

For the single-cell ATAC-seq data, Cell Ranger ATAC (v2.1.0, 10x Genomics, USA) was used to align the reads to the reference genome and calculate the binding activity across genomic regions. The binding sites and abundance information were read using Signac (v1.9.0, Stuart Lab, USA) (8), and the genomic regions were annotated accordingly. Nuclei with 3,000 – 100,000 unique fragments, transcription start site (TSS) enrichment score ≥8, and fraction of reads in peaks (FRiP) ≥0.20 were retained. Potential doublets were identified using the ArchR addDoubletScores() function with an expected rate of 5%, and nuclei with scores >1.0 were excluded prior to downstream analysis. The distribution of binding peaks across various genomic elements was computed and statistically analyzed. Data normalization was performed, and sample integration was achieved using the Harmony algorithm (v0.1.0, Korsunsky Lab, USA), followed by UMAP clustering. To integrate single-cell transcriptomic and ATAC-seq data, the FindTransferAnchors function was used to calculate transfer anchors, and single-cell ATAC-seq data were annotated based on the single-cell RNA-seq data. The R software was then used to calculate the distribution proportions of each cell type within the single-cell ATAC-seq data.

Using the Seurat package’s FindAllMarkers function, cell-type-specific regulatory peaks were identified (9). Overrepresented motifs were found by FindMotifs function. Subsequent functional enrichment analysis of these marker genes was conducted using the ClusterProfiler (v4.6.0, Yu Lab, China) and GSVA (v1.46.0 Hanzelmann Lab, Germany) packages to identify cell types closely related to the study’s focus (10). The distribution of these motifs across cell types and genes was analyzed using single-cell ATAC-seq data.

The scRNA-seq and scATAC-seq data were integrated using advanced computational approaches to dissect the cellular heterogeneity and regulatory programs in high-risk B-ALL. Differential gene expressions and chromatin accessibility analysis were performed to identify markers associated with treatment response and relapse. We constructed regulatory networks by integrating scATAC-seq data with scRNA-seq data, using SCENIC to infer potential regulatory interactions between transcription factors and their target genes (11).

For subgroup analysis of HSC/MPP cells, the relevant clusters were extracted from the integrated dataset. Cells were re-normalized using SCTransform, incorporating donor ID and percentage of mitochondrial gene expression as covariates. We selected the top 1,500 highly variable genes and conducted PCA, using the first 15 principal components for UMAP visualization and Louvain clustering. Batch effect correction was applied using the Harmony algorithm, with donor ID specified as the batch variable, and performed prior to dimensionality reduction. This approach ensured that downstream clustering reflected true biological heterogeneity within the HSC/MPP compartment.

We employed a robust bioinformatic pipeline that includes quality control, doublet detection, data normalization, dimensionality reduction using tools like PCA and t-SNE, and clustering algorithms to classify cells into distinct subtypes. Pathway enrichment analysis and gene set variation analysis (GSVA) were conducted to interpret the biological significance of the observed changes. Using the inferCNV algorithm (v1.0, Broad Institute, USA) for CNV analysis, with normal pro-B cells as references. “High CNV” refers to greater variability in chromosomal expression profiles as estimated by inferCNV, relative to reference Pro-B cells. Both gains and losses contribute to the CNV score.

All the visualizations, including UMAP plots, violin plots, heatmaps, dot plots, and bar charts, were generated using Seurat (v4.3.0), ggplot2 (v3.4.0), and ComplexHeatmap (v2.14.0) packages in R. Motif logos and accessibility tracks were plotted using ArchR and Signac.

Statistical analyses were performed in R unless otherwise stated. Comparisons between two groups were conducted using the Wilcoxon rank−sum test. For multiple comparisons, adjusted p−values were calculated using the Benjamini-Hochberg method to control the false discovery rate (FDR), with significance thresholds set at adjusted p < 0.05.

To construct the landscape of cell type-specific open chromatin features and gene expression profile at the single cell level, we performed sc-RNASeq and scATAC-seq analysis on PBMCs samples from six high risk pediatric B-ALL patients with and three healthy controls (CK1, CK2, CK3) to explore their cellular composition (Supplementary Figure S1). PBMCs were chosen due to their accessibility and reduced invasiveness, especially in pediatric patients undergoing frequent procedures. These high-risk patients had all undergone the standard CAML (Cyclophosphamide, Cytosine arabinoside, 6-Mercaptopurine, and Pegaspargase) chemotherapy regimen. Specifically, we selected three patients who achieved complete remission (CR1, CR2, CR3) and three who experienced non-remission (NCR1, NCR2, NCR3) following two courses of intensified treatment.

After an initial quality control assessment and doublet removal, we obtained a total of 76458 cells, with 38773 cells from healthy donors and 37685 cells from B-ALL patients, median number of estimated cells were 6861, and the ATAC Median high-quality fragments per cell were 5071. After performing normalization and log1p transformation, highly-variable gene selection, dimensionality reduction, batch correction, and Leiden clustering, cells originating from samples were separately annotated into distinct broad cell types and visualized via Uniform Manifold Approximation and Projection (UMAP). We identified 9 main clusters in parallel according to the gene profile. These clusters were annotated as CD8+ NKT-like cells, classical monocytes, HSC/MPP cells, memory CD4+ T cells, naïve B cells, naïve CD4+ T cells, naïve CD8+ T cells, pre-B cells, pro-B cells, and progenitor cells based on the expression of known markers for each cell type (Figures 1A, B) (12–15).

We then analyzed the relative proportions of cell clusters (Figures 1C, D). The cellular composition of the ALL remission group (CR group) is similar to that of the normal group (CK), with naïve CD4 T cells and CD8+ NKT-like cells being predominant. However, the cellular composition of the non-remission group (NCR group) is significantly different from the CR group, characterized by the highest proportion of Pro-B cells (p<0.0001), a marked increase in HSC/MPP cells (p=0.0001), and a significant decrease in naïve CD4 T cells (p=0.0002) and naïve B cells(Figure 1E).

The chromatin accessibility profile of individual cells is captured by scATAC-seq. Cell-type-specific chromatin accessibility profiles are relatively unknown. So, we used Seurat’s label transfer method to predict scATAC-seq cell types based on annotated scRNA-seq data. The detailed analysis methods were executed as previously described. Comparison between scATAC-seq cell-type predictions obtained by label transfer and curated annotations of unsupervised clusters indicates that all major cell types were present in both datasets and that scATAC-seq is comparable to scRNA-seq in the detection and assignment of cell identities (Figures 2A, B). Subsequently, we investigated the differences in chromatin accessibility among CK, CR and NCR groups (Figure 2C). In the CK group, naïve CD4 T cells were predominant, followed by classical monocytes and HSC/MPP cells (Figure 2D). In the CR group, CD8+ NKT-like cells were dominant, which is in line with the findings from scRNA-seq. Relative to the CK group, the CR group exhibited a notable increase in the proportion of memory CD4+ T cells and CD8+ NKT-like cells, alongside a significant reduction in the proportion of naïve CD4 T cells. The cellular composition of the NCR group diverged markedly from the other groups, with Pro-B cells predominating and accounting for more than half of the cells. In contrast to the CR group, the NCR group showcased a significant upsurge in the proportion of HSC/MPP cells (p<0.0001), naïve B cells(p<0.0001) and Pro-B cells (p<0.0001), coupled with a substantial decrease in the proportion of CD8+ NKT-like cells (p<0.0001), classical monocytes (p<0.0001), memory CD4+ T cells (p<0.0001), naïve CD4 T cells (p<0.0001), and naïve CD8 T cells (p<0.0001) (Figure 2E).

Given the most pronounced differences in the relative proportions of HSC/MPP cells and Pro-B cells between the CR group and the NCR group, and the concordance of findings from scRNA-seq and scATAC-seq multi-omics analyses, we concentrated our investigation on these two specific subgroups. We analyzed the differentially expressed genes (DEGs) in HSC/MPP cells and Pro-B cells between the NCR and CR groups. The gene ontology (GO) analysis of DEGs in HSC/MPP revealed that the B cell receptor signaling pathway and antigen receptor-mediated signaling pathway were upregulated in NCR (Figure 3A), while lymphocyte differentiation and T cell receptor signaling pathway were downregulated in NCR (Figure 3B). GO analysis of DEGs in Pro-B revealed that pathways governing the regulation of transcription, gene expression, and biosynthetic processes were enhanced in the NCR group (Figure 3C), contrasted with a suppression of signaling pathways involved in DNA replication within the NCR group (Figure 3D).

In order to further characterize the intrinsic structure and potential functional subtypes, we performed spectral clustering on HSCs cells using UMAP. As a result, we identified three stable subclusters in HSCs cells (Figure 4A). The CK group was mainly composed of subcluster 1, the CR group was predominantly made up of cells from subcluster 2, and the Pro-B cells of the NCR group were found in subcluster 0 (Figure 4B). Subcluster 0 included cells expressing TCF4, EBF1, ERG, AL589693.1, and CRIM1 genes (Figure 4C). Subcluster 1 included cells expressing genes like RPS27, RPL41, RPS15A, RPS18, and RPL13. Subcluster 2 included cells expressing genes such as RASA3, CD247, EML4, SKAP1, and ANKRD44. The TNF signaling pathway and T cell receptor signaling pathway were enriched in subcluster 2, while allograft rejection pathway and Notch signaling pathway were enriched in subcluster 0 (Figure 4D).

For the spectral clustering of Pro-B cells, five subclusters were identified (Figure 5A). The CK group was predominantly composed of subcluster 1, along with some cells from subcluster 4 (Figure 5B). The CR group was mainly made up of cells from subcluster 2, while the Pro-B cells from the NCR group belonged to subclusters 0 and 3. Each subcluster expressed its unique set of characteristic genes. The relative expression levels of cluster-specific markers in the five cellular subclusters of Pro-B have been presented in Figure 5C. RPS29, B2M, RPL41, and RPS21 were expressed across subclusters 0, 1, 2, and 3, but exhibited the highest expression in subcluster 1. In addition to RPS29, B2M, RPL41, and RPS21, subcluster 0 also included cells expressing the genes NEIL1, AC007384.1, and CRIM1. Subcluster 2 further included cells expressing genes such as HMGB2, PTPRC, SKAP1, and TUBA1B. Subcluster 4 was relatively distinct from the other subclusters, primarily consisting of cells expressing genes like GP1BB, NRGN, PF4, OAZ1, and PPBP. GSEA enrichment analysis showed that the B cell receptor signaling pathway was enriched in subcluster 0, while the cell cycle, DNA replication, and mismatch repair pathways were enriched in subclusters 3 and 2 (Figure 5D). Gene ontology enrichment analysis further revealed distinct functional programs among the HSC/MPP and Pro-B clusters (Figures 4E, 5E).

Copy number variation (CNV) analysis has been widely used in scRNA‐seq to investigate disease evolvement and development. To identify the cellular origin of malignant cells, we utilized the inferCNV algorithm to analyze the chromosomal CNV levels between different cell types (Figure 6A). As expected, the CNV levels in in the NCR group were significantly higher than those in the CR group (Figure 6A). Within the NCR group, HSS/MPP cells and pro-B cells exhibited higher CNV levels than other cell types.

We further analyzed the CNV levels in the three subclusters of HSC/MPP cells and the five subclusters of Pro-B cells. The results indicated that subcluster 0 of HSC/MPP cells had the highest CNV levels (Figure 6B), while in Pro-B cells, subgroup 0 exhibited the highest CNV levels, followed by subgroup 3 (Figure 6C).

Compared to other cell types, HSC/MPP cells in scATAC-seq had a large number of differential peaks of chromatin accessibility (Figure 7A). With transcription factor motif analysis, we found that the binding sites of FOSL2, FOS, FOSL1, BATF, JUN, JUNB had higher accessibility in the HSC/MPP than in other cell types (Figures 7B, C). However, only FOSL2 had a significant level of expression in the HSC/MPP group, while the expression levels of other genes were relatively low. When comparing HSC/MPP cells between NCR and CR groups, no differential peaks in binding site accessibility were observed (Figure 7B). Furthermore, expression patterns of representative immune-related genes across cell types were also assessed (Figure 7D).

The analysis of the Pro-B subgroup in scATAC-seq showed that the binding sites of EBF3, EBF1, MAZ, ZNF148, KLF15, NHLH1, and KLF5 had higher accessibility in Pro-B cells than in other cell types (Figures 8A-C). We further analyzed the expression differences of these genes in scRNA-seq and found that EBF1 and ZNF148 were highly expressed in Pro-B cells (Figure 8D). When comparing Pro-B cells between NCR and CR groups, the differential peaks in NCR were mainly downregulated, which means they were predominantly enriched in the CR group (Figure 8E). We found that the accessibility of KLF15, SP1, SP3, NRF1, ZBTB14, EGR1, EGR3 binding sites in Pro-B was higher in the CR group than in the NCR groups (Figure 8F). Among them, SP1, SP3, NRF1 also showed higher expression levels in scRNA-seq in the CR group compared to the NCR group, which is consistent with scATAC-seq (Figure 8G).