Single-cell Transcriptional Atlas of Human Hematopoiesis Reveals Genetic and Hierarchy-Based Determinants of Aberrant AML Differentiation Open Access

Supplementary Methods outlining computational analyses in this study

Supplementary Note 1: Comparison of transcriptional cell states and functionally-defined HSPC subsets Supplementary Note 2: Quantification of differentiation state abundance in bulk RNA-seq

Supplementary Tables S1 to S21, with Table Legends included

Supplementary Figure S1. Projection of sorted HSPC populations onto BoneMarrowMap A) Projection of purified multilineage (F1) and Mk/Ery restricted (F2, F3) sub-fractions of multipotent progenitor (MPP), common myeloid progenitor (CMP), and megakaryocyte-erythroid progenitors (MEP) profiled by bulk RNA-seq in Notta et al 2016 (44). B) Projection of scRNA-seq profiles purified with functionally informed gating of lymphoid-primed multipotent progenitor (LMPP), granulocyte-monocyte progenitor (GMP), and multi-lymphoid progenitor (MLP) populations from Karamitros et al 2018 (45). C) Projection of scRNA-seq profiles from hematopoietic stem cells (HSCs) purified using different markers of human HSCs. D) Proportion of transcriptional HSCs and MPPs projected from scRNA-seq data sorted with varying levels of immunophenotypic purity. Within each sample in this analysis, transcriptional HSCs and MPPs are defined using a stringent confidence threshold >0.8 based on K-nearest neighbours (KNN) classification after projection of query cells onto BoneMarrowMap. The combination of markers used, tissue source, and dataset name are also depicted. E-L) Projection of scRNA-seq profiles of immunophenotypically defined HSPC populations from Zhang et al 2022 (12), spanning E) HSC, F) MPP, G) LMPP, H) MLPs, I) Pre-B/NK, J) MEP, K) CMP, and L) GMP.

Supplementary Figure S2. Functionally distinct CD38 fractions within human GMPs and MEPs A-C) Re-analysis of index-sorted scRNA-seq profiles from Velten et al 2017 (8) based on CD38 surface expression within defined progenitor subsets. A) Gating of megakaryocyte-erythroid progenitor (MEP), common myeloid progenitor (CMP), and granulocyte-monocyte progenitors (GMP) populations within the CD34+CD38+CD10- compartment. B) Gating and scRNA-seq projection of CD38-mid and CD38-high GMP fractions. C) Gating and scRNA-seq projection of CD38-mid and CD38-high MEP fractions. D-F) Single-cell colony assays with index-sorting of CD38-mid and CD38-high fractions within defined E) GMP populations and F) MEP populations from umbilical cord blood. The index-sorting approach ensures that the recorded surface immunophenotype of each cell, including surface CD38 levels, can be linked to the functional colony output of that cell. G) Single-cell colony formation outcomes within CD38-mid and CD38-high GMPs and H) surface CD38 levels of individual GMPs that did and did not form colonies recorded by index-sorting. I) Single-cell colony formation outcomes within CD38-mid and CD38-high MEPs and J) surface CD38 levels of individual MEPs that did and did not form colonies recorded by index-sorting. P values from a two-tailed Wilcoxon rank-sum test are shown. Box plots indicate the range of the central 50% of the data with the central line marking the median. Whiskers extend from each both to 1.5x the interquartile range.

Supplementary Figure S3. Identification of gene expression programs through consensus non-negative matrix factorization (cNMF) A) Unsupervised discovery of 48 gene expression programs (GEP) across the human hematopoietic scRNA-seq atlas by cNMF. 36 GEPs correspond to hematopoietic cell states (cell state programs) and 12 GEPs correspond to other cellular processes (activity programs). The relative activity of each program within each hematopoietic cell state is depicted through a heatmap of AUC values. B) Transcription factor (TF) regulons associated with each GEP based on a composite score representing both the correlation of GEPs with each TF regulon score as well as the significance of overlap between genes defining each GEP and each TF regulon. C) Key pathways significantly associated with activity GEPs. D-G) Cell cycle plot of S vs G2M scores within hematopoietic stem and progenitor cells, depicting D) cell cycle assignments and E) S-phase, F) G2-phase, and G) M-phase program expression. H-I) UMAP of hematopoietic cell states depicting H) cell cycle assignments and I) S-phase, J) G2-phase, and K) M-phase program expression.

Supplementary Figure S4. Pseudotime analysis and characterization of erythroid and monocytic differentiation A-C) scRNA-seq projection and validation of hematopoietic pseudotime values in query data from Roy et al 2021 (10). A) Pseudotime values calculated de novo on hematopoietic stem and progenitor cell transcriptomes from Roy et al 2021 (10) using Monocle3. B) Pseudotime values inferred through BoneMarrowMap projection of HSPC transcriptomes from Roy et al 2021 (10). C) Comparison of de novo and projected pseudotime values across cell states annotated in Roy et al 2021 (10). Average values for each cell state are depicted and the size and color of each cell state correspond to the number of cells within that state. The linear regression line as well as r and P values from Pearson correlation are shown. D-F) Characterization of erythroid differentiation. D) Pseudotime score along erythroid differentiation. E) Transcription factor regulon activity along erythroid differentiation pseudotime. F) Gene expression changes along erythroid differentiation pseudotime. G-I) Characterization of monocyte differentiation. G) Pseudotime score along monocyte differentiation. H) Transcription factor regulon activity along monocyte differentiation pseudotime. I) Gene expression changes along monocyte differentiation pseudotime.

Supplementary Figure S5. Characterization of B-lymphoid and dendritic cell differentiation A-C) Characterization of B cell differentiation. A) Pseudotime score along B cell differentiation. B) Transcription factor regulon activity along B cell differentiation pseudotime. C) Gene expression changes along B cell differentiation pseudotime. D-F) Characterization of plasmacytoid dendritic cell (pDC) differentiation. D) Pseudotime score along pDC differentiation. E) Transcription factor regulon activity along pDC differentiation pseudotime. F) Gene expression changes along pDC differentiation pseudotime. G-I) Characterization of conventional dendritic cell (cDC) differentiation. G) Pseudotime score along cDC differentiation. H) Transcription factor regulon activity along cDC differentiation pseudotime. I) Gene expression changes along cDC differentiation pseudotime.

Supplementary Figure S6. Characterization of megakaryocyte and granulocyte differentiation A-C) Characterization of megakaryocyte differentiation. A) Pseudotime score along megakaryocyte differentiation. B) Transcription factor regulon activity along megakaryocyte differentiation pseudotime. C) Gene expression changes along megakaryocyte differentiation pseudotime. D-F) Characterization of eosinophil / basophil / mast cell (Eo/Baso/Mast) differentiation. D) Pseudotime score along Eo/Baso/Mast differentiation. E) Transcription factor regulon activity along Eo/Baso/Mast differentiation pseudotime. F) Gene expression changes along Eo/Baso/Mast differentiation pseudotime. G-I) Characterization of neutrophil differentiation. G) Pseudotime score along neutrophil differentiation. H) Transcription factor regulon activity along neutrophil differentiation pseudotime. I) Gene expression changes along neutrophil differentiation pseudotime

Supplementary Figure S7. Characterization of broad differentiation states implicated in AML A) Single-cell composition analysis of 318 patient samples, specifically showing the correlation of the centered log ratio (CLR)-normalized relative abundance of each projected cell state along the hematopoietic hierarchy, spanning 318 patient samples. B) Consolidation of precise cell states into broad differentiation states. Non-negative matrix factorization (NMF) was used to identify groups of correlated cell states, and the top-weighted cell states driving each NMF component were combined into broad differentiation states. Weights for each cell state underlying each NMF-defined differentiation state are shown in the heatmap. C-D) Enrichment of normal hematopoietic C) cell state programs and D) cell activity programs within each acute myeloid leukemia (AML) differentiation state, depicted through area under the receiver operating characteristic curve (AUC) values.

Supplementary Figure S8. Expression of therapeutic targets and LSC signatures by AML differentiation state A) Differential expression (DE) results from one-vs-all comparisons between differentiation states. DE was performed using pseudo-bulk profiles for each differentiation state from each of the 318 patient samples. Results for therapeutic targets in AML are shown. The log fold change (logFC) statistic for each gene is depicted through the size and color intensity of each dot, wherein red denotes higher expression while blue denotes lower expression. For visualization, logFC values were capped between -2 and 2. DE results with P < 0.05 are shown, those with FDR < 0.05 are starred. B-D) Leukemia stem cell (LSC) signatures across AML differentiation states. Cells belonging to each differentiation state from each patient were pooled into pseudo-bulk profiles and LSC signatures were scored by GSVA. Scores are shown for B) genes enriched in functional LSC fractions from Ng et al 2016 (57) and C) Quiescent LSPC signature from Zeng et al 2022 (30). P values from a two-tailed Wilcoxon rank-sum test are shown. Box plots indicate the range of the central 50% of the data with the central line marking the median. Whiskers extend from each both to 1.5x the interquartile range. D) Projection of pseudo-bulk profiles from scRNA-seq data classified based on prior AML cell type annotations from van Galen et al 2019 (15) and Zeng et al 2022 (30). Labeled cells from three datasets (van Galen et al 2019, Abbas et al 2021, Wang et al 2024) (15,19,27) belonging to each AML cell type (LSPC-Quiescent, LSPC-Primed, LSPC-Cycle, GMP-like, ProMono-like, Mono-like, and cDC-like) were collapsed into a centroid for each patient sample. The number of cells that comprise each centroid are depicted by the point size.

Supplementary Figure S9. Representative samples from each AML differentiation pattern A) For each AML differentiation pattern defined from cell state composition analysis, projection results for each of two additional representative scRNA-seq samples are depicted. Patient ID and originating dataset are labeled accordingly. B) Aggregate profiles of each AML differentiation pattern. For each differentiation pattern, cells were sampled evenly across constituent patient samples until a total of 3,000 cells was obtained and plotted for each differentiation pattern.

Supplementary Figure S10. Single-cell composition analysis within defined genomic subtypes of AML Single cell composition heatmaps are shown for AML patients belonging to specific genomic subtypes of AML, including: (A) KMT2A-rearranged AML (n = 33), (B) CBFB::MYH11 driven AML (n = 14), (C) RUNX1::RUNX1T1 driven AML (n = 18), (D) NPM1 mutated AML (n = 56), and (E) AML with TP53 mutations or Complex Karyotype (n = 34). The heatmap depicts the centered log ratio (CLR) normalized abundance of individual cell states within each patient sample, clipped at 0 and 4.5 for visualization. Each sample is annotated based on the study of origin, specific disease diagnosis, age group, and AML hierarchy classification. Cell state heterogeneity, quantified by the Shannon diversity index, is also shown for each sample. When present, recurrent genomic alterations are also annotated for each scRNA-seq patient sample.

Supplementary Figure S11. Inference of AML differentiation state abundance from bulk RNA-seq profiles A) Depiction of broad differentiation states in AML, each spanning multiple cell states. AML differentiation states were defined by single-cell composition analysis from Supplementary Figure S7B. Cell states which do not belong to a defined AML differentiation state (due to lack of correlation across AML samples) are shown in grey. B) Marker genes specific to each differentiation state in AML. Cells belonging to each differentiation state from each patient were pooled prior to differential expression analysis and heatmap visualization. C) LASSO regression to estimate AML differentiation state abundance in 302 pseudo-bulk samples wherein cells from individual patients were aggregated together. Marker genes for each differentiation state were used as features for LASSO regression, and cross-validation results of inferred abundances by LASSO in held-out patient subsets are depicted for each differentiation state. P values from a two-tailed Wilcoxon rank-sum test are shown. Box plots indicate the range of the central 50% of the data with the central line marking the median. Whiskers extend from each both to 1.5x the interquartile range. D) Number of positively and negatively weighted genes comprising each LASSO regression model for predicting abundance of each AML differentiation state in bulk RNA-seq data. E) Predicted and observed differentiation state abundance across 302 AML pseudo-bulk samples. For each association, the linear regression line as well as r and P values from Pearson correlation are shown.

Supplementary Figure S12. Clinical and genomic correlates of AML differentiation state abundance Associations between inferred abundance of AML differentiation states with genomic alterations across AML patient samples profiled by RNA-seq. The strength of each association, quantified as the test statistic from a generalized linear model adjusting for cohort as a covariate, is depicted through the size and color intensity of each dot, wherein green denotes higher abundance and purple denotes lower abundance. Only associations with unadjusted P < 0.05 by likelihood ratio test are shown, associations with FDR < 0.05 are starred. Associations are depicted for each of A) FAB morphology (n = 411 patients across TCGA, BeatAML1, and Leucegene) (59–61), B) WHO 2022 classification (n = 1166 across 5 cohorts re-analyzed by Severens et al 2024, consisting of TCGA, BeatAML2, Leucegene, TARGET, and LUMC) (31,58–62), C) ICC 2022 classification (n = 1192 patients across 5 cohorts re-analyzed by Severens et al 2024) (31,58–62), D) Hierarchy classification from Zeng et al 2022 (n = 864 patients spanning TCGA, BeatAML1, and Leucegene) (30,59–61), E) Gene expression clusters from Severens et al 2024 (n = 1224 patients across 5 cohorts re-analyzed by Severens et al 2024) (31,58–62), F) Integrated RNA + Hierarchy + Methylation based classification from Ma et al 2025 (n = 812 patients spanning TCGA, BeatAML1, and Leucegene) (59–61,65), and G) Disease history (n = 487 patients from BeatAML2). H-J) Associations for specific mutation patterns spanning 1224 patients across 5 cohorts re-analyzed by Severens et al 2024 (TCGA, BeatAML2, Leucegene, TARGET, LUMC) (31,58–62). Comparisons are depicted for H) CEBPA mutation type (n = 99) versus CEBPA wildtype, I) FLT3 mutation type (n = 340) versus FLT3 wildtype, J) NPM1 mutation combination (n = 249 with normal karyotype) versus NPM1 wildtype.

Supplementary Figure S13. Clinical relevance of lymphoid and erythroid lineage biases A) BoneMarrowMap projections of representative AML patient samples with mono-allelic mutations in the RUNX1 gene. B-D) Erythroid vs myeloid involvement within 136 AEL samples from Iacobucci et al 2019 (72), depicted as the difference between the inferred Early Erythroid abundance and GMP abundance in each sample. These are depicted in association with C) driver mutation status and D) clinical risk category as defined in the original study. P values from a two-tailed Wilcoxon rank-sum test are shown. Box plots indicate the range of the central 50% of the data with the central line marking the median. Whiskers extend from each both to 1.5x the interquartile range. E-F) Principal component analysis of AML differentiation states within 357 normal karyotype adult AML samples. Heatmap columns in (F) represent AML samples ordered by principal component 2 (PC2), representing the Lymphoid vs Erythroid axis. For each AML sample, mutation status as well as estimated relative abundance of each AML differentiation state are shown. To the left of the heatmap, a bar plot with the feature loadings for PC2 is shown, representing the specific coefficients assigned to each AML differentiation state for calculation of PC2. G) Volcano plot of Pearson correlations between the relative abundance of each AML differentiation state and ex vivo drug sensitivities (area under the dose response curve x -1) from the BeatAML2 (31) screen for 123 normal karyotype AML samples. Associations of differentiation states with higher drug sensitivities are colored in red and associations with lower drug sensitivities are colored in blue.

Supplementary Figure S14. Analysis of engraftment patterns among CD34 and CD38 sorted fractions in AML AML cells from 74 primary patient samples from Ng et al 2016 (57) were sorted into four fractions based on CD34 and CD38 expression and each fraction was transplanted into NSG mice at varying cell doses to screen for leukemia-initiating cell activity (L-IC). Injected right femur engraftment at 12-week post-transplant was evaluated by flow cytometry. Fractions that initiated CD45+CD33+ leukemic grafts were scored as L-IC+ (dark blue), those that did not initiate grafts were scored as L-IC (light grey), and those that initiated grafts with both CD33+ myeloid and CD19+ lymphoid cells were scored as multilineage (green) and presumed to be normal or pre-leukemic. Patient samples were grouped together based on engraftment patterns across their CD34 and CD38 fractions. For samples in which three or more CD34/CD38 sorted fractions were capable of engrafting NSG mice, the resulting CD34 and CD38 immunophenotypes of patient-derived xenografts (PDXs) from each injected fraction were inspected and samples were categorized based on whether immunophenotypic differences were observed between the PDXs generated from each injected fraction.

Supplementary Figure S15. Examples of AML patients with stable xenograft immunophenotypes produced across primary CD34/CD38 sorted fractions. Example of AML samples with stable xenograft immunophenotypes derived from multiple CD34 / CD38 sorted fractions. For each sample, the CD34 and CD38 immunophenotype is shown for their primary sample as well as for xenografts derived from each sorted fraction. For each sample, hierarchy classifications are also reported based on deconvolution of bulk RNA-seq profiles as per from Zeng et al 2022 (30).

Supplementary Figure S16. Sorting strategy for primary and xenografted AML samples A) Gating strategy for primary AML samples sorted into immature CD34+CD38- and mature CD34-CD38+ fractions, pertaining to Figure 7 and Supplementary Figures S17, S18, and S19. Live (PI-) single cells were gated on CD45+CD3- prior to sorting based on CD34 and CD38 status. Sorted immature CD34+CD38- and mature CD34-CD38+ fractions were subject to xenotransplantation as well as scRNA-seq. B) Gating strategy for AML xenograft samples profiled for scRNA-seq, pertaining to Figure 7 and Supplementary Figures S17, S18, and S19. Live (PI-) single cells were gated on CD45+CD3- prior to scRNA-seq. Detailed information on antibodies used for fluorescence-activated cell sorting (FACS) can be found in Supplementary Table S21.

Supplementary Figure S17. Cell state classification in AML pt.90240 primary and xenograft samples A) Unsupervised UMAP embedding and AML cell type assignments of scRNA-seq from primary samples of AML pt.90240. B) Unsupervised UMAP embedding and AML cell type assignments from patient-derived xenograft (PDX) samples of AML pt.90240. C) Density plot of unsorted bulk, CD34+CD38- sorted, and CD34-CD38+ sorted primary fractions from AML pt.90240. D) Density plot of PDX generated from CD34+CD38- primary cells and PDX generated from CD34-CD38+ primary cells from AML pt.90240. E) BoneMarrowMap projection of primary AML fractions from (C). F) BoneMarrowMap projection of xenografted AML cells from (D). 

Supplementary Figure S18. Distinct LSC-driven leukemia cell hierarchies co-existing within AML pt.90394 A-F) Experimental workflow for identifying co-existing LSC-driven hierarchies as in Main Figure 7, applied to pt.90394. A) Immunophenotype and scRNA-seq cell state composition of the primary AML sample. B) scRNA-seq projection results of the primary AML sample. C) Immunophenotype and scRNA-seq composition for the primary CD34+CD38- sorted fraction, representing 0.17% of primary cells. D) Immunophenotype and scRNA-seq composition for the PDX sample derived from injection of the CD34+CD38- sorted fraction. E) Immunophenotype and scRNA-seq composition for the primary CD34-CD38+ sorted fraction, representing 16% of primary cells. F) Immunophenotype and scRNA-seq composition for the PDX sample derived from injection of the CD34-CD38+ sorted fraction.

Supplementary Figure S19. Cell state classification in AML pt.90394 primary and xenograft samples A) Unsupervised UMAP embedding and AML cell type assignments of scRNA-seq from primary samples of AML pt.90394. B) Unsupervised UMAP embedding and AML cell type assignments from patient-derived xenograft (PDX) samples of AML pt.90394. C) Density plot of unsorted bulk, CD34+CD38- sorted, and CD34-CD38+ sorted primary fractions from AML pt.90394. D) Density plot of PDX generated from CD34+CD38- primary cells and PDX generated from CD34-CD38+ primary cells from AML pt.90394. E) BoneMarrowMap projection of primary AML fractions from (C). F) BoneMarrowMap projection of xenografted AML cells from (D).