Authors: Kyle Ferchen (1Division of Immunobiology, Cincinnati Children’s Hospital Medical Center, Cincinnati, OH, USA), Xuan Zhang (1Division of Immunobiology, Cincinnati Children’s Hospital Medical Center, Cincinnati, OH, USA), Kairavee Thakkar (2Division of Biomedical Informatics, Cincinnati Children’s Hospital Medical Center, Cincinnati, OH, USA; 3Department of Pharmacology and Systems Physiology, University of Cincinnati College of Medicine, Cincinnati, Ohio, USA), Guangyuan Li (2Division of Biomedical Informatics, Cincinnati Children’s Hospital Medical Center, Cincinnati, OH, USA), David Bernardicius (1Division of Immunobiology, Cincinnati Children’s Hospital Medical Center, Cincinnati, OH, USA), Sidharth Sen (2Division of Biomedical Informatics, Cincinnati Children’s Hospital Medical Center, Cincinnati, OH, USA), Priyanka Rawat (2Division of Biomedical Informatics, Cincinnati Children’s Hospital Medical Center, Cincinnati, OH, USA), Andre Olsson (1Division of Immunobiology, Cincinnati Children’s Hospital Medical Center, Cincinnati, OH, USA), Sierra N. Bennett (1Division of Immunobiology, Cincinnati Children’s Hospital Medical Center, Cincinnati, OH, USA), Crystal Potter (4Division of Immunology, Allergy and Rheumatology, Department of Internal Medicine, University of Cincinnati College of Medicine, Cincinnati, OH, USA), Fred D. Finkelman (4Division of Immunology, Allergy and Rheumatology, Department of Internal Medicine, University of Cincinnati College of Medicine, Cincinnati, OH, USA), Josh Croteau (5BioLegend Inc. (a Revvity company), San Diego, CA USA), Samantha Morris (6Division of Gastroenterology, Hepatology and Endoscopy; Division of Genetics, Brigham and Women’s Hospital; Department of Systems Biology, Harvard Medical School, Boston, MA, USA), Harinder Singh (7Departments of Immunology and Computational and Systems Biology, Center for Systems Immunology, University of Pittsburgh, Pittsburgh, PA, USA), Nathan Salomonis (2Division of Biomedical Informatics, Cincinnati Children’s Hospital Medical Center, Cincinnati, OH, USA; 8Department of Pediatrics, University of Cincinnati, Cincinnati, OH, USA), H. Leighton Grimes (1Division of Immunobiology, Cincinnati Children’s Hospital Medical Center, Cincinnati, OH, USA; 8Department of Pediatrics, University of Cincinnati, Cincinnati, OH, USA; 9Division of Experimental Hematology and Cancer Biology, Cincinnati Children’s Hospital Medical Center, Cincinnati, OH, USA)
Categories: Article
Source: Nature immunology
A unified multimodal single cell framework reveals a discrete state model of hematopoiesis in mice
Authors: Kyle Ferchen, Xuan Zhang, Kairavee Thakkar, Guangyuan Li, David Bernardicius, Sidharth Sen, Priyanka Rawat, Andre Olsson, Sierra N. Bennett, Crystal Potter, Fred D. Finkelman, Josh Croteau, Samantha Morris, Harinder Singh, Nathan Salomonis, H. Leighton Grimes
Abstract
Large-scale, unbiased single-cell genomics studies of complex developmental compartments, such as hematopoiesis, have inferred novel cell states and trajectories; however, further characterization has been hampered by difficulty isolating cells corresponding to discrete genomic states. To address this, we present a framework that integrates multimodal single-cell analyses (RNA, surface protein and chromatin) with high-dimensional flow cytometry and enables semiautomated enrichment and functional characterization of diverse cell states. Our approach combines transcription factor expression with chromatin activity to uncover hierarchical gene regulatory networks driving these states. We delineated and isolated rare bone marrow Lin^−^Sca^−^CD117^+^CD27^+^ multilineage cell states (‘MultiLin’), validated predicted lineage trajectories and mapped differentiation potentials. Additionally, we used transcription factor activity on chromatin to trace and isolate multilineage progenitors undergoing multipotent to oligopotent lineage restriction. In the proposed model of steady-state hematopoiesis, discrete states governed developmental trajectories. This framework provides a scalable solution for isolating and characterizing novel cell states across different biological systems.
A fundamental challenge in developmental biology is understanding the hierarchical states within stem cell and progenitor compartments, along with the gene regulatory networks (GRNs) that drive cell fate decisions.^1^ GRNs involve transcription factors that act on regulatory elements like enhancers and silencers, influencing cell-type-specific gene expression patterns.^2–4^ Single-cell genomics technologies like scRNA-seq and ATAC-seq^5,6^ provide high-resolution data on transcriptome and chromatin states. However, due to technical limitations, these techniques often rely on inferential rather than experimental evidence to understand developmental relationships and cellular potentials. Newer methods, such as CITE-seq^7^, combine transcriptome profiles and surface protein data, but do not permit the isolation of discrete populations because these methods are destructive. Integrating CITE-seq with high-dimensional flow cytometry (e.g., InfinityFlow^8,9^) provides a powerful approach to isolate and study cell populations manifesting discrete transcriptional and chromatin signatures, thus offering insights into the developmental potential of progenitors and their underlying GRNs.
In hematology, flow cytometry has long been used to identify and isolate progenitor populations based on surface markers. This method has helped define key progenitor types like multipotent progenitors (MPP),^10^ megakaryocyte erythroid progenitors (MEP), common myeloid progenitors (CMP), granulocyte monocyte progenitors (GMP)^11^ and common lymphoid progenitors (CLP)^12^, forming the basis of the classical model of hematopoiesis. However, flow cytometry is limited by the number of detectable markers, and recent findings show that these progenitor populations are more heterogeneous than previously thought.^3,4,13^ New techniques, including scRNA-seq, have revealed a more complex view of hematopoiesis, suggesting that lineage commitment is a continuum rather than a series of stable, defined states^14,15^. More recently, the punctuated continuum model re-introduced the classical concept of stability across developmental trajectories, with pools of cells that are variably lineage specified.^14,16,17^
Novel bioinformatics approaches and clonogenic assays have helped identify intermediate progenitor states, such as MultiLin cells, which exhibit mixed-lineage gene expression and reside within traditional CMP and GMP gates.^2–4^ These findings challenge classical progenitor definitions, highlighting the importance of transitional states in lineage specification. To better understand hematopoietic progenitors, we combined multiomic single-cell methods^7,18^ with high-dimensional flow-cytometric profiling^8,9^ and developed a unified computational framework (termed “ChromLinker”) that integrates data from different sources to derive GRNs that reflect developmental trajectories. Analyzing transcription factor activities on accessible chromatin regions, identified key surface proteins, such as CD55 and CD371, as markers of critical lineage transitions in MultiLin progenitors. This framework provides a comprehensive view of progenitor heterogeneity and reveals how transcriptional and chromatin profiles influence lineage commitment and orchestrate steady-state blood production; forming the basis of a model of steady-state hematopoiesis in which discrete states govern developmental trajectories. Interactive analysis tools and datasets for the MarrowAtlas are available at https://altanalyze.org/MarrowAtlas/.
Results
A multimodal single-cell atlas resolves multilineage states
To create a comprehensive CITE-seq^7^ atlas of murine bone marrow, we first identified optimal gating strategies to enrich for rare hematopoietic stem progenitor cell (HSPC) states (Fig. 1a). This includes broad CD117^+^ progenitors (Miltenyi autoMACS), Lin^−^Sca1^+^CD117^+^ stem cells and multipotent progenitors (hereafter HSC-MPP, excluding CD150^−^CD48^+^ MPP3-MPP4-gate cells) and Lin^−^CD117^+^CD127^+^ lymphoid progenitors (hereafter CD127^+^ progenitors, Extended Data Fig. 1a,b). Informed by retrospective and new scRNA-Seq captures, we developed a gating procedure for previously defined markers of MuliLin cells that depleted neutrophil and monocyte specified populations (Lin^−^CD117^+^CD34^+^CD115^−^Ly6C^−^ MultiLin gate)(Extended Data Fig. 1g,h)^2,3^.
Previously, we molecularly titrated CITE-seq antibodies with antibody-derived-oligonucleotide tags (ADT) on human bone marrow by diluting them five fold and analyzing ADT sequence abundance in scRNA-seq^19^. Initially, we manually-titrated 65 TotalSeq-A antibodies on CD117^+^ enriched mouse bone marrow cells^2,20^ and then captured the progenitor gates above using the TotalSeq-A Universal Cocktail (Universal Mix v1.0, n=195) (Supplementary Table 1, 2). As we previously showed in human^19^, manual titration improved detection of many ADTs over a universal mix (Extended Data Fig. 2a,b). As before, we performed individual 5-fold titrations of 195 antibodies (Extended Data Fig. 2d), prioritized for an optimized CITE-seq panel based on known discriminative markers of progenitor populations as well as dose-dependent signals of the antibodies used for detection (Fig. 1c). The revised titrated CITE-seq panel of 103 ADTs, demonstrated improved cell-state re-classification accuracy when used to profile 63,000 new cells in HSC-MPP, CD127^+^, MultiLin, CD117^+^ gates (online Methods, Extended Data Fig. 2a–c, 2e and Supplementary Table 1, 2). To annotate the captured populations, we compiled a database of reference cell states for HSC-MPP^21,22^, dendritic cell^23^, basophil-mast cell commitment^24^, early lymphoid^25^ and myeloid intermediate cells^2–4^. Then we integrated RNA and ADTs across the database of published cell populations, with those from unsupervised clustering of each modality (ICGS2^26^) using scTriangulate^27^ (Extended Data Fig. 2f). scTriangulate identified 87 stable discrete cell states, as evidenced by modality specific contribution scores and cluster confidence (Extended Data Fig. 2g–j). This multimodal integration defined clusters from a combination of prior-defined cell states (n=31) and unsupervised transcriptome (n=40) or ADT (n=16) clustering (Extended Data Fig. 2g–k). To supplement this atlas, we used a well-based scRNA-seq capture method (HIVE; Methods) to analyze CD117^+^ progenitors, and Lin^−^CD117^+^CD34^+^CD115^−^Ly6C^−^ MultiLin populations along with gates enriched for CD125^+^ eosinophil and FcER1a^+^ basophil-mast-cell progenitors (Extended Data Fig. 4a; Ly6C-1, and Ly6C^+^ Eo-Trajectory gates)(~15,000 cells)(Fig. 1d4). We annotated the cells through gene-set enrichment against markers from our database and from the literature ^28–34^. We noted several clusters unique to MultiLin (ML) gate which we denote as ML-1a, ML-1b, ML-2, and ML-3, along with those with transcriptomes enriched for cell cycle (ML-CC) or MDP gene expression (ML-MDP)(Fig. 1e). We also noted rare populations that that share qHSC and macrophage marker genes that selectively appeared in the HSC-MPP gate and occurred at the same frequency as macrophages (Extended Data Fig. 2l). These cells possessed unique transcriptional and cell surface markers relative to qHSC and macrophages (Extended Data Fig. 2l). Flow cytometry revealed CD193 and CD115 macrophage markers on sparse cells within the HSC-MPP gate (Extended Data Fig. 2m), consistent with reported naturally occurring HSC-macrophage interactions^35,36^. To explore gene expression and ADT in the 87 clusters or annotate submitted datasets, we developed a series of online analysis tools in MarrowAtlas (see Data Availability; https://altanalyze.org/MarrowAtlas/).
High-dimensional flow cytometry reconciles gating strategies
We generated a high-throughput flow cytometry atlas using the InfinityFlow protocol^9,37^ to match the CITE-seq atlas. We profiled 140 cell surface markers in bone marrow CD117^+^ cells from Irf8^GFP^ transgenic mice^38^. The backbone contained Irf8^GFP^ and 22 surface markers selected based on prior studies^3,4,11–13,20,34,39–48^ (Supplementary Table 3). Infinity markers include 95 antibodies from the titrated CITE-Seq panel (Supplementary Table 1), and 10 transcription factor transgenic reporters (e.g., Gata1, Gata2, Pu.1, Gfi1, Myc). The ‘Initial InfinityFlow Object’ marker panel was sequentially optimized to better separate known lineages (i.e., basophil, mast, eosinophil) to produce a ‘Curated InfinityFlow Object’ (Fig. 2a, Supplementary Table 3). Initial analysis with the software pyInfinityFlow^9^ yielded a UMAP embedding that resembled the structure found in the CITE-seq data (Fig. 2a).
To correspond the cells in our InfinityFlow atlas with prior studies, we performed in-silico gating using FlowJo to emulate defined hematopoietic progenitor gating schemes ^3,4,11–13,20,34,39–48^ (Fig. 2b). While gating long/short-term HSCs, MPP2, MPP3, and MPP4^46^ was largely consistent across studies, subdivisions within GMP and CMP (i.e., preGM, IG2) varied significantly across studies. For example, MPP1^48^ (which neighbors early megakaryocytic and erythroid specified cells)^13^, and the MPP4 compartment based on CD48 and CD62L expression^20^ were clearly delineated in the Infinity object. In contrast, In silico CMP and GMP gates^11^ evidenced substantial heterogeneity. Within CMP-GMP gating, CD150 identified megakaryocyte-erythroid biased cells^13^, transgenic fluorescent reporters for IRF8 and Gfi1 proteins identified the IG2^4^ population with monocytic and neutrophil potentials, and Ly6C expresion marked neutrophil or monocyte specified progenitors.^44,45^ We could also clearly detect CD125^+^ committed eosinophil progenitors (EoP)^42^, FCER1A^+^ basophil-mast cell progenitors (BMCP)^41^, and CD127^+^ common lymphoid progenitors (CLP)^12^. Thus, we integrated published progenitors within the InfinityFlow atlas, with the exception of CMP and GMP (which appear to be gates instead of populations).
Analyzing the InfinityFlow atlas, we noted that while CD34 persists on some myeloid progenitor states (e.g. eosinophil^42^ and neutrophil^3^ progenitors), Lin^−^Sca^−^CD117^+^ MultiLin cells were CD27^hi^ (Fig. 2c). Moreover, most bipotential lineage populations (except the ILC-NK progenitor branch) downregulated CD27 (Fig. 2c). Therefore, we used a Lin^−^Sca1^−^CD117^+^CD27^+^ gate to isolate MultiLin progenitors for the rest of the study. Thus, in silico gating of the InfinityFlow atlas enabled consolidation of a diverse set of historically delineated gating schemes for hematopoietic progenitors and established a foundation for integration with their CITE-seq and TEA-seq defined transcriptional and chromatin states. The Flow Cytometry Standard file for the final InfinityObject can be re-analyzed using standard commercially-available flow cytometry data analysis software (see Data Availability).
Integrated atlases permit isolation of cluster-defined cells
CITE-seq ADT values and InfinityFlow signals have distinct distributions and sensitivity. Assuming surface protein expression rankings remain consistent^49,50^, we tested whether the percentiles of surface protein expression remained consistent across platforms using a two-step, percentile-based normalization strategy (Fig. 3a). First, we mapped the CITE-seq antibody signal intensities for the gated HSC-MPP, MultiLin or CD127^+^ populations to the signal intensities of corresponding percentiles within the total CD117^+^ population. Next, we generated complementary spline functions to map InfinityFlow signals with their percentile ranks within the CD117^+^ compartment. This can be visualized and calculated as the area under the curve for the kernel density estimation (KDE) functions for the signal of each antibody. This two-step process enabled the inference of expected InfinityFlow fluorescence intensity for any given CITE-seq ADT value (Fig. 3a).
Using the software cellHarmony^51^ re-implemented in Python3, we transferred the transcriptome-defined CITE-seq cluster labels into the InfinityFlow atlas (protein features), resulting in well-separated populations (Fig. 3b). In addition to cell populations, capture gates from CITE-Seq generally correspond to the matching gates from InfinityFlow (Fig. 3c). Projection of the CITE-Seq directly within the InfinityFlow UMAP embedding enables direct correspondence between RNA expression and transcription factor reporter abundance, illustrating the fidelity of our alignment strategy (Fig. 3d, Supplementary Table 4). Comparison of the CITE-Seq projected labels with historical flow cytometry gating schemes by InfinityFlow (Fig. 2b) found a strong concordance for specified cell-states (i.e., erythroid, megakaryocytic, monocytic, dendritic, neutrophil), while MultiLin clusters overlapped with earlier delineations of CMP, preGM and GMP, suggesting our approach has the potential to resolve new cell states and gating solutions (Extended Data Fig. 3).
To rigorously test label transfer into the InfinityFlow object, we developed an in silico isolation strategy called flow-sort-scRNA-seq validation (“FS-scRNA-seq”) in which gates are predicted to enrich for target cell populations then utilized to capture populations for scRNA-Seq and mapped back to the atlas to evaluate enrichment. Given clear transcriptomic signatures^3,24,44^, we first isolated granulopoiesis states, and optimized surface marker combinations for sorting (Extended Data Fig. 4a). InfinityFlow-guided gates resulted in higher enrichment for the neutrophil subsets (Fig. 3e) than in silico gating (Extended Data Fig. 3). Within the Ly6C^−^Irf8^−^ GMP-gate fraction CD9 marked eosinophil and basophil-mast-cell lineage specification (Extended Data Fig. 4a). Eosinophil commitment exhibited a clear developmental progression as CD125^+^Ly6C^+^CD106^−^ EoP1 became CD125^+^Ly6C^+^CD106^+^ EoP2. While both EoP1 and EoP2 were absent in 10X Genomics captures, they were marked by a lineage-specific CRE reporter (Epx^Cre^ROSA^tdTomato+^)(Extended Data Fig. 4b–e) and were found in well-based captures albeit with low UMI counts (Extended Data Fig. 4f–i).
Next, we replicated gating schemes that sub-fractionated the CMP gate into progenitors specifying to granulocyte-monocyte or megakaryocytic-erythroid linages^13^. FS-scRNA-seq validation revealed that preMegE^13^ and erythroid/megakaryocyte gates^13^ were enriched for corresponding CITE-seq clusters (Fig. 3e). The preGM^13^ gate was highly heterogeneous and contained MultiLin clusters (Fig. 3e). To further enrich MultiLin clusters, we developed the algorithm Ab-MarkerFinder that iteratively optimizes surface marker combinations (based on unique expression patterns) for each target population (Methods). We applied this algorithm to all cell populations to predict optimal markers for isolation (Supplementary Table 5). These refined gates improved the isolation of several MultiLin populations (Extended Data Fig. 5a–5f). We achieved over 50% enrichment for the ML-1b cluster (Fig. 3e), while refined gates targeting specific MPP, MEP, BMCP and IG2 populations delivered high purity in flow-sorting-scRNA-seq validation (Fig. 3e). Thus, our InfinityFlow-derived strategy surpassed the resolution of historical gates to isolate unique populations spanning multipotent and oligopotent states.
To benchmark these correspondences, we took advantage of the fact that scRNA-seq captures from FS-scRNA-seq validated gates (Fig. 3e) could be used to provide labeled data in cases in which the transcriptome profile of specific flow cytometry gates is already known (Fig. 3f). By selecting populations with mutually exclusive transcriptome states (Extended Data Fig. 5g,h), we could quantify the success rate of different strategies to project labels from CITE-seq to InfinityFlow (Extended Data Fig. 5i). KDE mapping normalization plus cellHarmony label transfer outperformed the other tested techniques^49,52–54^ for transferring labels between CITE-seq and InfinityFlow (Fig. 3g). Thus, our framework for label transfer, which accounted for technological differences, outperformed alternative algorithms and signal normalization approaches.
Deep learning maps transcription factor activity
To delineate GRNs we profiled open chromatin regions within HSC-MPP- and MultiLin-gated cells using TEA-seq^18^, a trimodal assay that measures nascent transcriptomes, surface epitopes and chromatin accessibility (Fig. 4a and Extended Data Fig. 6). To define TEA-seq cluster analogs of the CITE-seq atlas, we applied the harmonypy^55^ workflow, and identified TEA-seq clusters corresponding to 57 of the 87 CITE-seq clusters with highly concordant marker genes (Extended Data Fig. 6c, Extended Data 7a–c). We validated the success of label transfer with Spearman correlations of marker gene rankings (Extended Data Fig. 7d, 7e). For the ATAC-seq cell states with sufficient read depth (32 of the 57) we inferred transcription factor activity using the ChromBPNet neural-network-modeling framework^56^. This enabled transformation of TEA-seq-cluster chromatin-accessibility profiles into base-pair-resolution ‘contribution scores’ that predict the importance of each DNA base pair to chromatin accessibility (measured as peak counts) within the given region (Fig. 5a, Extended Data Fig. 6). Short DNA sequences with similar contribution score profiles (seqlets) were clustered into contribution weight matrices (CWMs) using TF-MoDISco^57^ and assessed for their similarity to known transcription factor DNA-binding motifs using the CIS-BP2 database^58^. The CWMs were then scanned across all open chromatin regions to presumptively identify transcription factors that control accessibility and gene expression within each cell state. This analysis involved pairwise correlations of seqlet contribution score values for a given open-chromatin region and expression of a target gene located within the same topologically-associating domain (TAD)^59^. TF-gene interactions were prioritized based on the strength of Pearson correlation between seqlet contribution score and CITE-seq gene expression values across the 32 mapped TEA-seq clusters (Methods). Seqlets were found to cluster together with patterns specific to groups of CITE-seq clusters (e.g. Ets, RFX and C2H2-ZF factors in Sca^+^ HSC-MPP clusters, bZIP and Ets factors in MultiLin, neutrophil and monocytic clusters, and GATA and bHLH factors in erythroid and megakaryocytic clusters; Fig. 4b,c).
To define cell-state-specific regulatory interactions within the stem/progenitor compartment, we restricted ChromBPNet predicted transcription factors to targets overlapping top cluster marker genes. For this analysis, we used the CITE-seq RNA expression data paired with the seqlet contribution scores of linked open chromatin regions and known transcription factor DNA-binding motif annotations to prioritize transcription factor-target gene links (Methods, Extended Data Fig. 6; Supplementary Table 5). Cytoscape network visualization of the predicted TF-gene interactions in various cell states was consistent with known developmental functions of the transcription factor Hlf in stem cells, Gata2, Gata1, Myc and Nrf1 in megakaryocytic-erythroid progenitors and Spi1 and Irf8 in monocytic-dendritic cell progenitors (Fig. 4d). Cell-state-specific models highlighted the known activity of Hlf in qHSCs, which lacked lineage-specific transcription factor activities (Fig. 4e, Extended Data Fig. 8a). MultiLin states ML-1b and ML-2 showed step wise priming of Pu.1 and Irf8 activity, which was more pronounced in IG2-proNeu1 and ML-MDP (Fig. 4e, Extended Data Fig. 8a). In contrast, ML-3 showed nascent Gata1 and Gata2 activity, which was retained in BMCP and became more pronounced in MEP (Fig. 4e, Extended Data Fig. 8a). Thus, cluster-specific GRNs were constructed from base-pair resolution contribution scores, the expression of transcription factors and their target genes.
To generate a comprehensive map of transcription factor activities across the 32 identified TEA-seq clusters, we utilized a composite activity Z-score by integrating transcription factor expression, target gene expression, and the regulatory contribution of the transcription factor to its putative target genes in each of the 32 clusters (Methods). This analysis revealed a distinctive set of transcription factors that were more active in HSC-MPP and MultiLin cells (Fig. 4f). Importantly the regulatory network inferred activities of SPI1, GATA2, GATA1, IRF8, and MYC corresponded with their InfinityFlow transcription factor-reporter expression values across the 32 cell states (Fig. 4g–j and Extended Data Fig. 8b). Thus, activity scores were well supported by their InfinityFlow reporter expression values after integration by cluster label propagation (Fig. 3d,e).
Surface proteins predict transcription factor activities
We sought to identify surface proteins that reported the activity of transcription factors underlying hierarchical lineage specification. When we aggregated activity scores for individual genes encoding dynamically expressed surface proteins in the MultiLin clusters, we noted that CD55 was expressed on early basophil, mast cell, erythroid and megakaryocytic progenitors and that Cd55 expression was linked to the activity of Gata1 and /or Gata2 transcription factors, while CD371 was expressed on early eosinophil, neutrophil, monocyte and dendritic cell progenitors and Clec12a expression was linked to the activity of Spi1, Irf8, Cebpa and Cebpe (Fig. 5a, Extended Data Fig. 9b). Using our CITE-Seq expression data, CellOracle^60^ in silico knock-out (KO) of Gata1 or Gata2 downregulated Cd55 expression and shifted trajectory away from CD55^+^ ML-3 (Fig. 5b–c, Extended Data Fig. 9b, Supplementary Table 6). Conversely, KO of Spi1, Irf8, Cebpa, or Cebpe downregulated Clec12a expression and shifted trajectory away from CD371^+^ ML-2 (Fig. 5b–c, Extended Data Fig. 9b, Supplementary Table 6). Thus, while inducing the expression of CD55 or CD371, these transcription factors act dynamically within discrete MultiLin cells to regulate the specification of erythroid or myeloid cell fates.
To detemine whether the genomic regions implicated in the assembled GRN (Fig. 4b–d) were unique or known, we tested their enrichment in published epigenomic profiling experiments.^61^ The majority (~68%) of seqlet positions in the GRN overlapped with previously identified candidate cis-regulatory elements (cCRE)(ENCODEv4^61^), across all 32 TEA-seq clusters (Extended Data Fig. 8c–d). To validate the identity of the transcription factors inferred to act on these cCRE (Extended Data Fig. 8e), we used GIGGLE indexing and search^62^ to query a collection of over 3,500 ChIP-seq datasets covering 347 transcription factors (CistromeDB^63^). We found a strong (~90%) enrichment between the transcription factor families corresponding to seqlets (Gata1, Gata2, Irf8, Spi1, Cebpa, and Cebpe) and their ChIP-seq binding sites (Extended Data Fig. 8f and 8g). Restricted by TADs^59^, analysis of cCREs identified candidate binding positions for CEBP factors (Cebpa and Cebpe), as well as ETS-IRF composite elements^64^ (Irf8 and Spi1) correlated to expression of Clec12a (Fig. 5d, 5e), while binding sites for GATA factors (Gata1 and Gata2) correlated to expression of Cd55 (Fig. 5f). Thus, the transcription factor inferences were validated by published ChIP-seq data^63^, and their chromatin activity linked to expression of target genes.
Transcription factors induce hierarchical lineage restriction
Lineage priming was suggested to be dictated by the onset of multilineage transcription factor activity that is selectively reinforced and resolved in distinct patterns that determine particular lineages^4^. To determine the extent of mixed-lineage priming in each cell population, we computed the likelihood that a given cell state coincidently expresses defining markers for all committed lineage cell states (termed sc-Hrödinger). HSC had the lowest sc-Hrödinger scores, while MultiLin clusters had among the highest scores across CITE-Seq cell states (Extended Data Fig. 9a). Marker genes for neutrophils (e.g., Elane, Ptrn3, Ms4a3), eosinophils (e.g., Hdc, Cebpe, Ldhc), monocytes (e.g., F13a1, Ly6c2, Slpi) and dendritic cells (e.g., Ccr2, Rab7b, Ms4a4c) were primed in ML-2 cells, whereas marker genes for MkP (e.g., Pf4, Slc14a1, Gp5), erythroblasts (e.g., Klf1, Aqp1, Tspo2), basophils (e.g., Ms4a2, Alox5, Cd200r3) and mast cells (e.g., Scin, Gpr183, Kcnc1) were primed in ML-3 cells (Extended Data Fig. 9b). An alternative approach to identify fate transitions, Capybara^65^, identified lineage priming predictions concordant with sc-Hrödinger (Extended Data Fig. 9c).
Next, to infer lineage trajectories in our CITE-Seq atlas between predicted cell states, we applied the partition-based graph abstraction software (PAGA^66^). PAGA predicted 95 cell state interactions and prior-validated precursor-progeny relationships among bi-potential intermediate cell-states (i.e., IG2 to neutrophil and monocyte) and HSC-MPP to megakaryocytes (Extended Data Fig. 9d), suggesting that many inferred relationships are real^67^. To experimentally test the PAGA predictions, we used CellTag^68^, a lentiviral library with expressed barcodes and GFP, captured in single-cell transcriptomes to facilitate lineage-tracing and enable clonal tracking (Extended Data Fig. 9e). Both Sca^+^ (Lin^−^Kit^+^Sca^+^CD27^+^) and Sca1^−^ (Lin^−^Kit^+^Sca^−^CD27^+^; ‘MultiLin’) populations were sorted from C57/Bl6 bone marrow, transduced with CellTag vectors and cultured as techical replicates with BMEC-Akt^69^ immortalized endothelial cells for 48h. GFP^+^ cells were sorted and captured for scRNA-seq (Extended Data Fig. 9e) and the transcriptomes were annotated based on the CITE-seq atlas clusters using cellHarmony. Using a prior described clonal coupling score^70^, considering statistically significant interactions with evidence from both technical replicates (Methods), we identified 128 predicted cell state interactions, representing potential lineage trajectories (Fig. 6a, Online Methods). ML-1b had connections to itself and to ML-2, few ML-1b cells were identified (Supplementary Table 7 and Source Data). Similarly, few mast cells, EoP and basophil cluster cells were identified (Fig. 6b), indicating a bias in the developmental output or capture efficiency of specific populations in this in vitro system. CellTag interactions confirmed 63 of the 95 relationships predicted by PAGA (Extended Data Fig. 9d). When projected into the CITE-seq atlas UMAP (Fig. 6b), the CellTag interactions confirmed known lineage relationships and nominated new ones, as noted below (Supplementary Table 7). Newly-detected relationships that involved MultiLin states had a lower clonal coupling score (6) compared to relationships between known populations (15.6) (Supplementary Table 7), consistent with the transitory nature of the MultiLin states. Strong links (clonal coupling score >6) were observed between MkP, ERP, and ML-3 populations (Fig. 6b), and between MDPs, monocytes and neutrophil populations (Extended Data Fig. 9b), suggesting a divergence between CD55^+^ ML-3 and CD371^+^ ML-2 MultiLin states. Distinct clonal relationships were observed for ML-2 (BMCP, IG2-MP, MDP-2) and ML-3 (MkP-HSC, ERP-HSC)(Fig. 6b), suggesting they were distinct cell states, rather than arbitrarily delineated subsets. These clonal relationships were also independently observed in a prior generated inducible in vivo lineage-recording mouse dataset of hematopoiesis (DARLIN^71^)(Fig. 6b).
To address the functional potential of MultiLin populations, we used lineage-specific EpoR^Cre^ ROSA^−LSL-tdTomato^ (erythroid), Epx^Cre^ ROSA^−LSL-tdTomato^ (eosinophil), Mrp8^Cre^ ROSA^*-LSL-tdTomato^ (neutrophil) and CD11c^Cre^ ROSA^−LSL-tdTomato^ (dendritic cell) reporter mice. Introducing bone marrow from the latter mice into the InfinityFlow object (Fig. 6c). Notably, Cre activity was not detect in MultiLin cells (Fig. 6c). Individual cells sorted from CD55^+/−^CD371^+/−^ MultiLin gates (Extended Data Fig. 10a) were cultured for 5 days in StemSpan SFEM with cytokines (Methods) and assessed in terms of clonogenicity and expression of tdTomato (Extended Data Fig. 9f). Cells sorted from the CD55^−^CD371^−^ MultiLin gate gave rise to tdTomato^+^ cells marked by all four CRE (Fig. 6d), indicating they represented an early progenitor. Cells sorted from the CD55^+^CD371^−^ MultiLin gate, but not the CD55^−^CD371^+^ MultiLin gate gave rise to EpoR^CRE^Rosa^tdTomato+^ cells; whereas cells sorted from the CD55^−^CD371^+^ MultiLin gate gave rise to tdTomato^+^ cells marked by Epx^-*CRE^, Mrp8^−CRE^ and CD11c^−CRE^ (Fig. 6d). These data indicated that nascent and distinctive transcription programming within the diverse Lin^−^Sca^−^CD117^+^CD27^+^ MultiLin compartment selectively restricted their developmental potentials.
MultiLin states respond to perturbation
We developed a 27-color flow cytometry panel and visualized populations in the embedding space based on InfinityFlow in-silico gating (Fig. 7a). We found general concordance between the expected cluster content of in-silico-sorted InfinityFlow populations and their flow-sort-sRNA-seq analysis (Extended Data Fig. 10b). The lineage markers CD117, Sca1, CD27 and Ly6C enriched the MultiLin fraction (Fig. 7a); CD371 and Irf8 identified the initial split to eosinophil, neutrophil, monocyte and dendritic cell outputs (Fig. 7a); while CD55 revealed basophil, mast, erythroid and megakaryocyte cell outputs (Fig. 7a). Within the CD371^+^ populations, Irf8^lo^ cells could be divided based on Ly6C (neutrophils) or CD125 (eosinophils) expression (Fig. 7a) while CD115 or CD135 expression marked monocytes or MDP, respectively in the Irf8^hi^ fraction (Fig. 7a). The CD55^+^ cells could be split by differential expression of CD150, ITGB7 and CD41 for Meg/eryth lineage populatoins, while differential expressin of CD131, FCER1A and CD117 could distinguish basophil/mast lineage populations (Fig. 7a). The remaining CD55+CD150−CD131− gate contained an ITGB7+ population, which could be an erythroid-mast cell progenitor (EMaP), because they localized between erythroid and mast cell outputs in the InfinityFlow object (Fig. 7a) and ITGB7− population (Fig. 7a). Overall, this combination of markers provided a discrete separation of MultiLin fractions.
To examine the developmental regulators underlying the separable MultiLin states, we analyzed their transcription factor expression profiles using in-silico gating of the InfinityFlow object (Fig. 7b–d). The CD371^−^CD55^−^ ML1 gate cells expressed a mixture of antagonistic transcription factors, including intermediate levels of PU.1, IRF8, GFI1, GATA2 and GATA1 (Fig. 7e). This gating scheme detected previously-identified bistable switches^72,73^, including divergent PU.1 and GATA1 expression between ‘ML1’ to “ML2” and “ML3” gate cells, an IRF8+GFI1+ population neighboring IRF8^lo^ (eosinophil and neutrophil) and IRF8^hi^ (monocyte and DC) populations, and a GATA1^+^GATA2^+^ population neighboring GATA2^hi^ (basophil and mast cell) and GATA1^hi^ (megakaryocyte and erythroid) populations (Fig. 7e). To trace the development of the gated populations, we used marrow from the lineage-specific-CRE ROSA^−LSL-tdTomato^ reporter mice (as in Fig. 6), and index sorted single cells across the gates to assess CRE marking. The “ML1” gate cells gave rise to progeny marked by all four CRE reporters (Fig. 7f). EpoR^Cre^ ROSA^tdTomato+^ (erythroid) fate was enriched in the “MEP” and “ML3” gates (Fig. 7f); Epx (eosinophil) fate was enriched in “ML2” and “EoNP” gates (Fig. 7f); Mrp8 (neutrophil) fate was enriched in “ML2”, “EoNP”, “MP” and “MDP” gates (Fig. 7f); CD11c (DC) fate was enriched in “ML2”, “MP” and “MDP” gates (Fig. 7f), indicating the gating scheme resolved the divergent lineage potentials of the MuliLin cells by integrating cell surface markers, chromatin and gene expression states as well as expression of developmentally important transcription factors.
To determine whether MultiLin states respond to infection or stress, we analyzed their alterations in N. brasiliensis infection, which is known to induce a Th2-driven response, with the accumulation of basophils, mast cells and eosinophils in lungs and gut^74^. In Balb/c mice infected with Nippostrongylus brasiliensis (subcutaneous injecitoin of 750 L3 larvae as published^75,76^) full-spectrum flow cytometry and scRNA-seq capture of bone marrow cells prior to infection (D0) and at days 5, 7 and 10 after infection(Extended Data Fig. 10c). In comparison to D0 we detected a specific increase in the bone marrow CD55^+^CD371^−^ ML3 gate at day 5 and day 7 and a simultaneous decrease in the CD371^+^CD55^−^ ML2 gate (Fig. 7g,h), which coincided with an increase in bone marrow CD131^+^ basophil progenitors and FcER1a^+^ mast cell progenitors at day 7 (Fig. 7g,h). Differential gene expression analysis revealed a specific upregulation of genes in ML-3, BMCP, mast cell and basophil cluster cells starting at day 5 versus day 0 (Fig. 7i). In contrast, EoP cluster cells exhibited an upregulation of similar genes starting at day 3 (Fig. 7i), while the ML2 cluster cells (the putative EoP precursor) were reduced at day 3, and did not show the same transcriptional signal (Fig. 7i). The IL-5 receptor (CD125) was detected as specific marker of EoP (Source Data), consistent with the fact that infection-driven IL-5 production^77^ amplifies CD125^+^ EoP^78^, indicating that, while infection-induced IL-5 directly signaled to nascent EoP, the expansion of the ML3 gate cells that proceeded the increase in BMCP gate cells confirmed its role as a precursor to the basophil/mast cell fate. Thus, the specific responses of MultiLin populations to N. brasiliensis infection supported their developmental and physiological importance in the hematopoietic system.
To provide additional evidence for the role of ML-1, ML-2 and ML-3 cells across diverse in vivo perturbations (including genetic, cancer, infection and aging), we performed a retrospective analysis of 13 published scRNA-seq datasets^79–91^ spanning 22 distinct perturbations (see Methods). This analysis identified reproducible concordant shifts in ML-1 that corresponded to increased lymphoid production, in ML-2, corresponding with monocytic, dendritic and neutrophil production and in ML-3 corresponding with BMCP and ErP/MkP enriched outputs (Extended Fig. 10d–e). Thus, MultiLin subsets behaved as key determinants of biased lineage outputs of the hematopoietic system in the context of diverse perturbations spanning, infection, cancer and disease.
Discussion
In this study, we developed a unified framework that integrates multimodal single-cell analyses (CITE-seq, TEA-seq, and InfinityFlow) to isolate and characterize distinct hematopoietic progenitor populations. By combining transcriptional, chromatin, and surface protein features, we identified novel multilineage populations and mapped their functional potentials, offering a broadly applicable strategy for complex tissue systems.
This framework bridges classical surface-marker-based flow cytometry with genomic approaches, enabling the isolation and validation of progenitor states based on dynamically regulated transcription factors and stable molecular features such as open chromatin and surface proteins. We identified multipotent progenitors beyond traditional gates and characterized MultiLin cells based on their distinct gene regulatory networks (GRNs). Using this framework, we were able to isolate cells with the expected gene expression patterns, track the developmental potentials of different progenitor states and explore the dynamic GRNs underlying lineage specification. A public resource for exploration of tehese GRNs is provided (https://altanalyze.org/MarrowAtlas/).
One of the key observations from this study was the dynamic instability of MultiLin cells, which express co-regulated alternative-lineage determinants.^4^ This instability suggests a stochastic process of lineage restriction, as supported by prior findings in HSPCs^4,92^ and enteroendocrine progenitors.^93^ The MultiLin GRNs reflect dynamic activities of CEBP, IRF, and GATA family transcription factors, which cooperate or antagonize each other to regulate open chromatin regions.^94–96^ The resulting bifurcation induces the emergence of oligopotent progenitors, such as ML-2 and ML-3, which can be distinguished based on surface markers like CD371 and CD55. These results challenge previous models suggesting lineage bias within HSC-MPP and instead highlight the role of nascent transcription factor activity in MultiLin to promote lineage restriction.
We argue that hematopoietic progenitors exist in discrete, hierarchically organized states (multipotent, oligopotent, bipotent, and lineage-restricted) defined by their integrated transcriptional, chromatin, and proteomic features. The hierarchically nested GRN are the fundamental basis of discreteness and stability of developmental states within the hematopoietic compartment and lineage trajectories. The GRN enable dynamic control of hematopoietic outputs during homeostasis and in response to various stresses. The distinct GRN clearly illustrate the complex developmentally-diverse progenitor content of classical hematopoietic flow gates^10–12^, and argue against lineage commitment as a continuum.^14,15^ We do not fully exclude the possibility of a continuum of lineage progression with some stable states, as proposed by the punctuated hematopoiesis model.^14,16,17^ However, as opposed to pools of variably lineage specified progenitors^14,16,17^, we argue that the intermediate states identified here (e.g. MultiLin cells) represent more stable hierarchical developmental stages with distinct developmental potentials determined by their underlying GRN. These stable states were validated using cell barcoding and lineage tracing methods, providing strong evidence for their role in controlling developmental outcomes. Our data also suggest that clonal memory and gene regulatory mechanisms likely influence the production of different lineages, which could provide insight into clonal heterogeneity and lineage restriction in hematopoiesis suggested by the “parallel hematopoiesis” model.^97^ Indeed, single HSC transplant results in heterogeneous production, some of which is unilineage.^98^
Online Methods
Mice
All procedures were performed according to an Cincinnati Children’s Medical Center Institutional Animal Care and Use Committee approved protocol ( IACUC2023–0009).
Mice were housed on a 14hour light/10 hour dark cycle at 72^o^F, in ventilated micoisolator cages that are autoclaved sterile with irradiated isopad (medical grade cotton) bedding. Filtered reverse osmosis water pouches were made fresh daily. The diet was irradiated Teklad 2919. The facility is AAALAC accredited, and managed by specialized and ACLAM board-certified veterinarians.
Mice of both sexes between the ages of 4 weeks to 8 weeks were used for these studies. Some mice were obtained from the Jackson Labs (Bar Harbor ME); including C57/Bl6 (Strain #:000664), Balb/c (Strain #:000651), and transgenic C57/Bl6 mice including Irf8-GFP transgenic mice^38^(Strain #:027084), Plzf-GFPcre^99^ (Strain #:024529), Csf1r-EGFP^100^ (Strain #:018549), GFP-c-Myc KI^101^ (Strain #:021935), Ncr1gfp^102^ (Strain #:022739), Tcf7GFP flox^103^ (Strain #:030909), Ki67-RFP^104^ (Strain #:029802), ROSA-LSL-tdTomato^105^ (Strain #:007914), Epor-CRE^106^ (Strain #:035702), Mrp8-CRE^107^ (Strain #:021614), CD11c-CRE^108^ (Strain #:008068). Other transgenic mice were gifts of the Pu.1-YFP^109^, Gata1-mCherry^110^, Gata2-Venus^72^, Gfi1-tdTomato^111^, Tox-GFP^112^, Epx-CRE^113^, Ms4a3-CRE^45^, FUCCI2^114^.
Bone Marrow Cell Isolation and CD117 Enrichment
Mouse bone marrow cells were isolated from iliac crest, femur, and tibia bones of 6- to 8-week-old C57/BL6 male mice or Balb/cJ for N. brasiliensis infection modeling. Bones were crushed to obtain single cell suspension and CD117 enriched (130–097-146, Miltenyi Biotec) on Miltenyi AutoMACS (program Possel) according to manufacturer’s protocol.
Single-cell CITE-seq generation
CITE-seq was performed with 3 different ADT panel configuration of BioLegend TotalSeq-A (1) a customized 65-plex cocktail (manually titrated using flow cytometry), (2) a prototype 200-plex cocktail for antibody titration (Universal 1.0 BioLegend Catalogue Number: 199901), and (3) a final titrated and lyophilized 110-plex cocktail (Supplementary Table 1: available as BioLegend custom Cat: 900003823) with fresh CD135 and CD127 added as spike-in ADT.
The 65-plex cocktail was comprised of known hematopoietic stem and progenitor markers as an initial test of the technology. ADT conjugated antibodies were incubated with AF647-fluorochrome-conjugated oligo-dT oligonucleotide (IDT). Using CD117-enriched mouse bone marrow as input, concentrations were estimated from flow cytometric titration of the AF647 signal. After review of these data, final concentrations were recommended by BioLegend. All antibodies used and associated adjustment from vendor concentrations are indicated in Supplementary Table 1.
The 200-plex cocktail was provided as a lyophilized mix, uniform for each TotalSeq antibody, which was rehydrated at 8x concentration and subsequently diluted to five 4x, 2x, 1x, 0.5x, and 0.25x. 1x concentration was determined based on the performance of ADTs in the 65-plex panel or according to manufacturer’s recommendation. Five sorted populations were profiled with the 200-plex panel to comprehensively represent rare cell populations in the total bone marrow cells, CD117-enriched bone marrow cells, FS HSC-MPP cells (Lin^−^Sca1^+^CD117^+^ and exclusive of CD48^+^CD150^−^), FS MultiLin cells (Lin^−^CD117^+^CD34^+^CD16/32^+/low^CD115^−^CD11b^−^Ly6C^−^CD150^−^), and FS CD127^+^ (Lymphoid) cells. A unique HTO antibody was added to denote each titration concentration. Same cell population stained with different titration concentrations were pooled together after washing separately.
For the titrated 110-plex panel, a customized lyophilized cocktail was created by BioLegend using selected concentrations from the 200-plex titration experiment. The titrated 110-plex cocktail (with fresh CD135 and CD127 added as spike-in ADT) was applied to the same five populations except for total bone marrow cells.
In addition, each sorted cell population from various gates were stained with a unique HTO, washed and pooled. In each experiment, cells were washed on Laminar Wash MINI or HT2000 System (Curiox Biosystems Inc) with the following 25 cycles, flow rate of 10uL/s, initial volume of 55uL. Cells were captured using the Chromium X with Chromium Next GEM Single Cell 3ʹ kit v3.1 chemistry (PN-1000268, 10X Genomics).
Fluidigm scRNA-seq
Flow cytometry sorted cells purified from murine bone marrow were captured using the C1 Single-Cell Auto Prep System (Fluidigm) and underwent scRNA-seq analysis as previously described.^3^
TEA-seq
TEA-seq was performed on HSC-MPP gate and MultiLin gate cells following the published protocol^18^ with minor (1) using a digitonin concentration of 0.025% based on optimization results, (2) the addition of Protector RNAse Inhibitor (Roche 3335402001) in perm, wash and tagmentation buffers, (3) two ports were loaded to created duplicate technical replicates, and (4) 80ul instead of 160ul supernatant from cDNA amplifications were used to purify pre-amplified ADT. Briefly, 100,000 to 250,000 leftover cells from CITE-seq experiment were permeablized in digitonin perm buffer, washed and then processed according to the manufacturer’s protocol for Chromium Next GEM Single Cell Multiome ATAC + Gene Expression (PN-1000285, 10X Genomics).
Single cell library prep and sequencing
Library prep was performed according to protocols with minor modifications such as reduced PCR cycles for ADT and HTO library amplifications. We found pre-PCR ADT/HTO in our libraries were abundant enough that we generally reduced 2 cycles (ADT: 10–12 cycles, HTO: 9–10 cycles). Final transcriptome, ADT and HTO libraries were quantified and analyzed by Qubit dsDNA HS (1000268, Invitrogen), High Sensitivity DNA kit (5067–4626, Agilent Technologies) on 2100 Bioanalyzer (G2939BA, Agilent Technologies) and KAPA HiFi library quantification kit (KK4824, Roche). Dual-indexed transcriptome libraries were pooled and sequenced across multiple Illumina Novaseq 6000 S4 or X plus flow cells (Illumina) with PE150+10+10 or PE100+10+10 settings. Single indexed ADT and HTO libraries were pooled and sequenced on Illumina S2 flow cells with PE50+8 setting or with transcriptome libraries. BCL files were demultiplexed into fastq files for Cell Ranger input. “AT” was added to the end of RPI-x ADT i7 index (6bp) to match D70X_long HTO index (8bp). HTO and ADT FASTQ files were supplied as 3P feature barcode together with transcriptome FASTQ files into Cell Ranger V6.1.2 count pipeline. Transcriptome was mapped to mm10-v2 reference genomes for downstream analysis and visualization.
HTO calling and quality control
Cells were multiplexed using HTOs to distinguish CITE-seq ADT concentration or cell population. HTO barcode count matrices were obtained through the multimodal analysis workflow in Cell Ranger, prior to normalization (counts per ten thousand (CPTT)). Cell barcodes with > 30% of normalized reads assigned to multiple HTOs were annotated as doublets, with confident singlet assigned to cells with >40% of normalized reads assigned to a single HTO (HTO processing module of AltAnalyze v.2.1.4). Cells were further filtered based on the 7 mouse/rat isotype control antibodies counts (Source code) and performed quality control filtering in Seurat V4 by nFeature_RNA>500 & nCount_RNA >1000 & percent.mt<25. This QC step filtered 393,748 cells to 315,792 high-quality single cells in the initial titration dataset and 90,889 to 72,198 cells in the final titrated CITE-seq dataset. Individual library metrics are provided for all single-cell genomic experiments in Supplementary Table 2.
Multimodal analysis
All Cell Ranger produced count matrices underwent ambient RNA exclusion using the software SoupX^115^ with a contamination fraction of 15% and quality control filtering by HTO. Ambient corrected transcriptome counts and associated ADTs counts were supplied as input to the software TotalVI^116^ to obtain normalized and denoised ADT counts. To derive clusters from the initial titration CITE-seq datasets, the software cellHarmony was used to transfer labels from CPTT normalized expression centroids computed in author-provided labels from three prior published reference bone marrow atlases. Marker heatmaps were obtained using MarkerFinder in AltAnalyze using either single cells or combined donor pseudobulks for each scTriangulate cell population. For differential ADT analyses, we applied an empirical Bayes moderated t-test (FDR corrected). To quantitatively assess performance of the final titrated ADT mix relative to the original TotalSeq-A Universal 1.0 mix, we generated ADT mean expression values for all 87 scTriangulate clusters for all ADTs (TotalVI normalized) in both datasets and re-assigned cell identify to all cell-barcode using these centroids with cellHarmony (AltAnalyze 2.1.4). Performance was evaluated using Adjusted Rand Score function in the python scikit-learn library.
HIVE scRNA-seq
The Honeycomb HIVE CLX capture system (Honeycomb Biotechnologies, Inc) was applied using the manufacturer’s instructions to the following CD117-enriched bone marrow, sorted cells from basophil/mast enriched gates, and those from eosinophil enriched gates. Namely, CD117+ progenitors, and Lin-CD117+CD34+CD115−Ly6C− MultiLin populations along with eosinophil- and basophil-mast-cell-progenitor enriched gates (Extended Data Fig. 4a; Ly6C-1, Ly6C+ Eo Trajectory gates). Briefly, 30,000 cells of each population were loaded per HIVE by centrifugation and processed immediately without freezing. Libraries were quality checked and then pooled for sequencing on a Novaseq 6000 SP flow cells for PE50 using customized sequencing and indexing primers. Cells were called based UMI per cell knee plot and gene count matrix was generated by aligning to mm10-v2 for downstream analysis.
Decision-level cluster integration
The software scTriangulate (version 0.13.0) was applied to the optimized titrated CITE-seq compendium to define high-confidence multimodal single-cell clusters. As input for scTriangulate, clusters were derived from three separate sources as outlined 1) literature-centric supervised analysis, 2) independent transcriptome unsupervised clustering and 3) ADT unsupervised clustering. scTriangulate applies coalitional iteration with diverse single-modality or multi-modal input clustering solutions, to assess aggregate stability (Shapley value) of each overlapping annotation at a single-cell level. scTriangulate was run with 5 stability metric 1) Re-assign score, 2) Single Cell Clustering Assessment Framework (SCCAF) reclassification score and 3) Term Frequency-Inverse Document Frequency (TF-IDF) ranking the top first, fifth and tenth markers as separate scores (GitHub Repo). Among the 90 resulting scTriangulate clusters, three clusters labels and cells were excluded due to high mitochondrial gene expression. The final clusters were refined by computing centroids for these 87 clusters for 100 random cells per cluster, reclassifying all parental cells using the software cellHarmony in AltAnalyze v.2.1.4 (correlation cutoff > 0.3). Cell population names were initially derived from Fisher exact test enrichment compared to the literature-defined cell cluster labels and manual curation based on identified marker genes. Prior to scTriangulate, literature-centric cell population annotations were derived from the aggregate of prior murine progenitor single-cell cluster annotations, using reference-based label projection. Reference centroids or defining gene-sets were selected from studies in which reliable evidence exists to demonstrate progenitor lineage potential (CFU, index sorting, secondary transplant). Specifically, HSC-MPP subsets, dendritic cell, basophil/mast cell commitment, early lymphoid, and myeloid intermediate centroids were derived from each study^2–4,11,21–25^ and used to identify the best matching cells in our pre-titration scRNA-seq compendium using the software cellHarmony (top 50 scoring cells per centroid). The top-selected unique cells from all references were combined to produce an aggregate reference, using defined marker genes (MarkerFinder, AltAnalyze). The combined multi-study reference centroids were applied back to the optimized titrated CITE-seq dataset, using cellHarmony centroid classification (default options). To derive unsupervised RNA and ADT clusters in this dataset, we applied ICGS2 clustering (default options) in two rounds, to define subclusters for the principle identified lineages in AltAnalyze (Ensembl version 72 BioMart database). For each modality, sub-clustering was applied to initial ICGS2 clusters binned according to flow-gate inferred lineages (HSPC, MultiLin, CD117, CD127 enriched). From 28 initial transcriptome defined clusters, 58 ICGS2 sub-clusters were identified. From 13 initial ADT defined clusters, 93 ICGS subclusters were identified. Full source annotations are provided in Synapse (syn60529836).
Cell fate prediction analysis
sc-Hrödinger is a python module in the software AltAnalyze (https://github.com/nsalomonis/altanalyze/blob/master/stats_scripts/multiLineageScore.py) developed to compute the probability that a given cell state is consistent with mixed-lineage priming (metastability) based on the coincidence expression of lineage-defining markers from presumed uni-potential committed progenitor cell states. sc-Hrödinger quantifies the degree to which a cell simultaneously expresses markers for multiple lineages (in >1 defined reference cell states, in >25% of cells of that cell state). Multipotent cell states are frequently defined by priming to one more lineages. For markers per reference cell states, for a given cell, a sc-Hrödinger is computed as the mean of the gene-level calls for expression all lineage marker genes (1 or 0, if CPTT > 0 for each gene). Unique non-overlapping marker genes are defined using the AltAnalyze MarkerFinder algorithm. In this analysis, MarkerFinder was run on 11 selected uni-potential committed progenitor cell states representing the major cell lineages identified in our CITE-seq analysis (aHSC, MkP, mast cell progenitor, basophil progenitor, eosinophil progenitor, monocyte progenitor, early neutrophil progenitor, early lymphoid, early erythroid, early conventional dendritic cell, early plasmacytoid dendritic). The top 100 marker genes per representative cell states were provided (Pearson Rho ranked) along with an expression file filtered to these genes (CPTT CITE-seq normalized counts) and associated cellbarcode to cluster annotations.
Capybara Cell Identity Analysis
Single-cell RNA sequencing data from the titrated CITE-seq experiment was processed and log-transformed from CPTT values after soupx ambient RNA correction. Cell identity was quantified using Capybara (v1.0). First, a reference dataset was constructed using specific end-state clusters. Second, cell identity was measured across all of the captured cell states. Quadratic programming scores representing the continuous cell identity assignments were extracted from the query dataset to generate a heatmap visualization of z-score normalized QP values. Cells exhibiting characteristics of multiple end states were curated using the Capybara function multi.id.curate.qp and the frequency of “Multi-ID” states for each input CITE-seq cluster identity were plotted as a percentage.
MultiLin Cell State Correlations with Perturbation
To determine the functional outputs associated with perturbations in distinct MultiLin cell states, we performed a retrospective analysis of 13 published scRNA-seq datasets^79–91^ spanning 22 distinct perturbations pre-leukemic models, infection, inhibitors, knockouts, aging and disease) in single-cell bone marrow RNA-seq (GSE227026, GSE248396, GSE228562, GSE191147, GSE252833, GSE223632, GSE197407, GSE209742, GSE147729, GSE189217, GSE264087, GSE236407, GSE235798). Label transfer to scTriangulate CITE-Seq centroids (syn66721893) was performed using an optimized version of cellHarmony for direct use with 10x Genomics Chromium h5 and mtx format files (https://github.com/SalomonisLab/altanalyze3). Significant associations between the cell frequency of distinct MultiLin subsets and downstream differentiation outcomes were quantified using standardized chi-square residuals [(observed – expected)/sqrt(expected)] on the cell frequency counts table (R chisq.test function). The rank correlation of the standardized residuals were calculated to determine the extent to which there was an association between over(under)-representation of pairs of cell types. Multiple testing and dependencies among the p-values was accounted for using the Benjamini-Yekutieli (BY, R p.adjust function) FDR procedure, where a BYfdr ≤ 0.05 was considered statistically significant. This analysis yielded MultiLin associations with 64 down-stream populations that were visualized as a heatmap.
RShiny app development
Two RShiny apps were developed from the final Titrated ADT CITE-seq (1) an Azimuth reference atlas for label transfer of uploaded scRNA-seq datasets and (2) a ShinyCell interactive browser for exploration of single-cell populations, RNA and ADT expression and distinct dataset covariates. These apps are available at https://altanalyze.org/MarrowAtlas/. The multimodal Azimuth bone marrow reference RShiny interface was built following the Azimuth v0.4.6 instructions (https://github.com/satijalab/azimuth) using the neighbors from the titrated RNA data restricted to the top MarkerFinder marker genes (syn66721893). CITE-seq RNA counts were scaled and normalized as CPTT with clusters defined for four different annotation levels based on the multimodal scTriangulate clusters. A ShinyCell viewer for the healthy bone marrow compendium was generated using a formatted h5ad counts matrix with corresponding sample/cell-level metadata for both the titrated dataset (syn66721894). Associated RShiny and h5ad creation scripts are linked to the datasets in Synapse.
Infinity Flow
Infinity Flow captures were performed using a 5-laser spectral cytometer (Aurora System, Cytek Biosciences). An initial backbone of 22 selected pan-lineage cell-surface protein markers were used for all InfinityFlow captures to enable regression to impute 95 additional fluorochrome-conjugated antibody signals, and 18 transgenic fluorescent reporter signals. All Flow Cytometry fluorescence intensities were normalized using the logicle transformation as implemented in the pyInfinityFlow package using default parameters. Regression was carried out using pyInfinityFlow to apply the XGBoost algorithm^54^, enabling the imputation of 113 features flow cytometry features, overlapping with the optimized CITE-seq panel, including the original backbone. As ground truth for XGBoost, we leveraged Flow Cytometry of ~100,000 cells stained with one of the Infinity Markers and all of the backbone antibodies. We trained marker specific XGBoost models for the measured Infinity Marker on 80% of cells and tested the accuracy of the regression using mean squared error between the predicted Infinity Marker expression and the ground truth Infinity Marker expression on the remaining 20% of cells “held out” from training, using the default settings of pyInfinityFlow. For the prediction of surface proteins, CD117-enriched bone marrow cells were sampled from mice with the Irf8-GFP transgenic reporter to improve prediction accuracy. The target sample on which the final regression models were applied included staining for CD131-PE to not impute this signal and carry no error from prediction. This yielded an initial InfinityFlow object based on the 22-color backbone. To identify poor performing imputed markers to optimize the original 22 marker backbone, we assessed imputation variability (mean squared error) or selected markers that were highly expressed in cell populations with poor separation in the obtained UMAP embedding. This analysis nominated 27 additional markers able to distinguish multiple rare cell populations, including early specifying basophil, mast, and eosinophil cells. This updated panel was applied to a new target sample from Irf8-GFP CD117-enriched bone marrow samples. Overlapping markers between the initial InfinityFlow object and the new target population were used to generate regression models to predict all remaining signals from the initial InfinityFlow object onto the curated panel capture. The resulting curated InfinityFlow object was used for all downstream processing. In-silico gating was done by importing the InfinityFlow object into FlowJo as an FCS file. FCS files are available in Synapse (syn60529836).
Supervised KDE Mapping
To normalize the CITE-seq ADT signals to best approximate the signals of analogous surface proteins observed in flow cytometry, reference spline functions were built to map the percentile of expression level for each feature to that feature’s intensity in the InfinityFlow space. This creates a map between rank and signal intensity that is assumed to be the same for any sample drawn from the same population (CD117-enriched murine bone marrow cells). The signal intensity is assumed to be conserved between the different populations captured by CITE-seq after TotalVI batch correction (CD117-enriched, HSC-MPP gated cells, MultiLin gated cells, and CD127+ gated cells), so another spline function is created to map the signal intensity of the CD117-enriched population ADT intensity for each feature to its percentile rank among cells. These CITE-seq specific spline functions can then be applied to each cell and for each feature in all sorted populations to map their signal intensity to the predicted percentile across the distribution observed for CD117-enriched bone marrow. That percentile then serves as input to the analogous InfinityFlow reference spline function to map to the InfinityFlow signal intensity. This two-step, percentile-based normalization strategy is conceptually similar to landmark registration, which aligns datasets based on corresponding features.^117^
FS-scRNAseq Capture
Surface marker expression as well as CITE-seq atlas cluster label transfer to the InfinityFlow object were used to nominate new gating schemes. Initial gates targeted specified progenitor states for a first round of FS-scRNA-seq capture. These gates were applied to a Sony MA900 sorter. In a second round of FS-scRNA-seq capture, gates from Pronk et al. were applied using a Cytek Aurora CS. Since the Sony MA900 and Aurora CS operate on different bit systems, optimization measures were taken when setting up the CS to ensure accuracy (a comparison sort was run in parallel on Aurora CS and BD FACS Aria). A third round of FS-scRNA-seq capture used a new brute force algorithm called Ab-MarkerFinder, which iteratively applied the MarkerFinder algorithm (pyInfinityFlow) and in silico gating to nominate an optimal gating strategy (details below)(Supplementary Table 5). These were applied in silico with visualization over the 2-dimensional UMAP space to curate markers for the final gating strategy. All final gating strategies were integrated into a single panel and applied on the Cytek Aurora CS. We subsequently performed FS enrichment of all cell populations using conventional flow cytometry. The majority of these populations were captured using the Chromium X with Chromium Next GEM Single Cell 3ʹ kit v3.1 chemistry (PN-1000268, 10X Genomics), along with HSC-MPP and CD127+ produced CITE-seq captures.
Ab-MarkerFinder
To nominate novel gating schemes for populations of interest from the curated InfinityFlow object, we developed a new brute-force algorithm that could be trained with millions of cells. Briefly, MarkerFinder is applied for the prospective population of interest to identify the top 5 correlated and top 5 anti-correlated features. For each of these features, the signal is split over n (default n=100) evenly spaced amounts after logicle normalization of fluorescence intensity. For each of these values, a candidate gate in both positive and negative direction is tested for improvements to both purity and yield of the prospective population. The best gate (determined by improvement of purity), picked among the 10 tested markers, is then used to filter the data. The process is then re-iterated with the filtered data for the prospective population. The iteration continues until either the desired purity (default 90%) or a yield limit (default 50%) is reached.
Benchmarking of ADT to InfinityFlow label transfer approaches
To assess the accuracy of KDE mapping relative to prior established multimodal label transfer approaches, prior described and novel hematopoietic gating strategies with corresponding single-cell captures were used at a benchmark dataset. This benchmarking dataset was restricted to gates that were observed to be mutually exclusive when visually assessed using the combined CD117-enriched bone marrow InfinityFlow object. Following the same QC filtering protocol and SoupX correction approach as the CITE-seq data, the top 60 unique marker genes per each of the 23 FS isolated cell populations using the software MarkerFinder (250 random cells per cluster). To obtain ADT profiles for each cell population, corresponding cells in the final titrated CITE-seq data identified using cellHarmony transcriptome mapping with two separate parameters to ensure rigor (centroid- and community-based alignment).
Using these matched InfinityFlow and CITE-seq ground-state predictions, we evaluated mapping of the Flow-defined CITE-seq populations to the corresponding InfinityFlow in silico flow gate defined populations. These analyses consider the 113 matching CITE-seq ADTs (or corresponding gene mRNA) and InfinityFlow cell-surface markers. To assess the precision of KDE mapping with cellHarmony alignment relative to potential orthogonal strategies, we tested 4 distinct multimodal label transfer approaches (Seurat Bridge-integration, MARIO-CCA, Harmony-KNN, XGBoost) ^52–55^ and CyCombine^49^ feature normalization. Bridge integration was implemented with the source FS-scRNA-seq sorted populations and associated labels (input), the titrated CITE-seq (bridge) and the filtered InfinityFlow object (target). Harmony-KNN was applied using harmonypy Python library to the combined CITE-seq and InfinityFlow cell surface marker intensity matrix, without KDE mapping, followed by a k-nearest neighbors classification model built using scikit-learn. A python XGBoost classification model was fitted to the 23 labels in the CITE-seq matrix and queried against the InfinityFlow matrix. As ground truth for XGBoost, we leveraged Flow Cytometry of ~100,000 cells stained with one of the Infinity Markers and all of the backbone antibodies. We trained the XGBoost model for the measured Infinity Marker on 80% of cells and tested the accuracy of the regression using mean squared error between the predicted Infinity Marker expression and the ground truth Infinity Marker expression on the remaining 20% of cells “held out” from training, using the default settings of pyInfinityFlow. For MARIO-CCA, the MARIO algorithm was similarly applied with canonical correlation analysis (CCA) vectors. For the cellHarmony using community-alignment, was applied with the pyInfinityFlow associated library using either KDE mapping, no mapping or CyCombine normalization. The ARI score was used to assess accuracy of the obtained labels from each separate approach to the ground-state in silico flow InfinityFlow population annotations.
Integration of CITE-seq and InfinityFlow and UMAP Projection
After normalizing the ADT intensity values with the supervised KDE mapping approach, cellHarmony (pyInfinityFlow package) was applied to transfer the scTriangulate defined labels to the InfinityFlow dataset using community-alignment. Upon transferring these labels, the InfinityFlow dataset was subsampled to a size of 1e5 cells, distributing the sampled cells as evenly as possible among the cluster labels. This subsampled InfinityFlow set was used to fit a 2-dimensional UMAP projection using the 27 base fluorescence signals of the captured InfinityFlow object. This UMAP was then applied to the entire InfinityFlow dataset using the transform function to apply find embedding coordinates for all cells. The resulting coordinates and cluster labels were saved to the InfinityFlow object FCS file. To co-embed the CITE-seq data, a k-nearest-neighbors (KNN) model was used to identify the top 4 nearest neighbors of the CITE-seq cells to the 1e5 sampled InfinityFlow cells, and the UMAP position was taken to be the centroid position of those 4 nearest neighbors in the embedded space (Consensus KNN). To assess the correlation between CITE-seq mRNA count and InfinityFlow transcription factor reporter fluorescence intensity, the Pearson correlation was calculated between the CITE-seq mRNA count for those transcription factors against the centroid fluorescence intensity of the InfinityFlow reporter expression value.
In-vitro HSPC Culturing Conditions
Bone marrow HSPCs were cultured using StemSpan SFEM (StemCell Technologies), enriched with 50ng/mL SCF, 10ng/mL IL-3, 10ng/mL IL-6, 2U/mL EPO, 50ng/mL IL-11, 10ng/mL IL-5, 50ng/mL TPO, 15ng/mL G-CSF, 15ng/mL GM-CSF, 10ng/mL IL-7 (PeproTech, Rocky Hill, NJ USA).
CellTag Lineage Mapping
pSMAL-CellTag-multi-V1 plasmid DNA library was validated by Sanger sequencing.^118^ The lentiviral vector used in these studies was manufactured by the Vector Production Facility at Cincinnati Children’s Hospital Medical Center, as described briefly below. The target vector, Delta 8.9, VSVG, and pRSV-Rev plasmids were packaged in DMEM using PEI transfection reagent. Transfection reagent was removed after 6 hours. Viral supernatant was harvested 44- and 68-hours post transfection, clarified with a 0.45 um filter, and processed with Gamma Gold Clarification, XT5 Scale Ion Exchange Chromatography, and Tangential Flow Filtration to remove contaminants from cell culture.
Virus was titrated on Lin^−^Sca^−^CD117^+^CD27^+^ MultiLin cells and titer was chosen on transduction efficiency of 80%. CellTagging was performed separately on two cultures of Sca^+^ (Lin^−^Sca^+^CD117^+^CD27+ gate) and two cultures of Sca1^−^ (Lin^−^Sca^−^CD117^+^CD27+ gate ‘MultiLin’) marrow progenitors. These technical replicates were processed independently, including single cell captures and library generation. 50,000 (Sca+ or Sca−) cells per well were plated at 50uL in a 96-well plate, then spinfected (10,000 RPM, 45min at room temperature), and then incubated overnight in a 37C incubator with 5% CO2. Cells were recovered the next day and transferred to a 24-well plate with bone marrow epithelial cells (BMEC-Akt^69^) to potentially preserve self-renewal capacity. BMEC-Akt were plated at 80% confluency 48 hours prior to co-culture. 0h starts when transduced HSC/progenitor cells enter co-culture with BMEC-Akt. After 48 hours of co-culture, transduced cells were sorted for GFP^+^CD45.2^+^Ly6C^−/low^ to both gate away BMEC-Akt1 cells, and avoid over-representation of neutrophil output cells. Technical replicates were processed separately to create four 10X 3’ scRNA-seq libraries (Sca+1, Sca+2, Sca−1, Sca−2).
CellTag Clonal relationship inference
The CellTag^118^ associated scRNA-seq was aligned to our CITE-seq compendium cell atlas annotations using the software cellHarmony as described above. UMAP coordinates were projected from the CellTag into the CITE-Seq embeddings using a custom python script (Synapse: syn53237568). To reconstruct lineage relationships from the CellTag, we assessed the cell identities shared in each clone and calculated a comprehensive clonal coupling score between two cell types, as described below.
Clonal coupling score.
We calculate a clonal coupling score^119^ between two cell types. We first calculate 1D45Cij the number of observed unique clonal barcodes, or clones, shared between two cell types, i and j. We calculate an expected number of shared clones between the two cell types using the oij values in the Oij matrix, where i,j∈c=[c1,...cn], and c=[c1,...cn] is the set of cell types in the following Eij=(ΣcOci)(ΣcOjc)/(Σc,cOc,c)
Clonal coupling score between two cell types is calculated as the ratio of observed matrix Oij over the expected matrix Eij.^119^ The clonal coupling score, thus, is a symmetric adjacency matrix which can be visualized as a network. Firstly, the clonal coupling scores were centralized. Lineage relationships for which the clonal coupling score was or greater than 30 was set to a maximum value of 30. Similarly, lineage relationships for which the clonal coupling score was or less than 1/30 was set to a minimum value of 1/30. Lastly, the clonal coupling scores were log2-normalized. Note that in contrast to the approach by Weinreb et al^119^, clonal coupling scores were not calculated as the median observed/expected matrix ratio from randomized trials of sampled cell-type specific cells. We specifically avoided the latter to account for smaller clones from rarer cell populations.
Clonal relationship inference.
To leverage the confidence acquired through technical replicates, we consider all the lineage relationships between two cell types that satisfy one of the two 1) log-normalized clonal coupling score > 0 in both replicates, and 2) number of shared clones > 4 in both replicates. This criterion assumes that a lineage relationship is unlikely to occur out of random chance if it is observed in both the replicates, thereby increasing the specificity of a lineage relationship. Additionally, the criterion accounts for relationships that are penalized by clonal coupling score for having the number of observed barcodes too close to the expected number of barcodes.
Lineage Trajectory Inference
Log-scaled normalized gene expression (CPTT) of only the MarkerFinder genes (n = 2,893) from the cells from the CITE-seq was used to create an anndata object using the scanpy package. Two different approaches, diffusion pseudotime (DPT) and PAGA, from the scanpy package were used to infer pseudotemporal ordering of the hematopoietic stem cells. Highly variable genes were determined with the default parameters. UMAP coordinates of titrated CITE-seq data were evaluated in PAGA. For PAGA, cell populations were restricted to those quantified using CellTag to enable direct comparison.
Diffusion Pseudotime (DPT).
Neighborhood graph of cells was computed on 30 principal components (‘n_pcs’ parameter) with the number of neighbors (‘n_neighbors’ parameter) set to 30, as recommended by the scanpy authors.^120^ Diffusion map was computed with the default settings. Cell populations were restricted to those detected by CellTag to draw direct inferences. To calculate the diffusion pseudotime, a cell annotated as “MPP1-G1” was set as the root cell for pseudotime (most inferred primitive HSC-MPP cluster captured by CellTag). In our analysis, the first MPP1-G1 barcode in the list (“AGGGAGTAGCTGCCTG-1.AS_CITE_HSC”) was set as the root cell. Diffusion pseudotime was calculated using the function scanpy.tl.diffmap and the pseudotime values for each cell are projected on the UMAP.
PAGA.
PAGA^66^ graph was computed using the connectivities from the neighborhood graph (described above) and annotated single-cell clusters as nodes. “MPP1-G1” was set as the root node (‘root_key’ parameter). DPT values are used to guide the PAGA graph (‘use_time_prior’ parameter). Finally, only the edges with weight greater than 0.6 in the PAGA are visualized on the UMAP.
Network similarity metrics
Jaccard Index.
The Jaccard index is scaled from 0 to 1, where 0 means the least similarity in the edges of two networks. It is the ratio of intersection of edges between two networks divided by the union of edges between the two networks.
Hamming–Ipsen–Mikhailov (HIM) metric.
HIM metric is a previously published network similarity score^121^ to assess the similarity between the topologies of the two networks, taking into account differences in edge weights and degree (the number of edges leaving a node or cell type) distributions. As indicated by the equation below, HIM distance is a linear combination of the normalized Hamming distance, which accounts for differences in edge weights, and the normalized Ipsen–Mikhailov distance, which accounts for the similarity in degree distributions. The Ipsen–Mikhailov has parameter γ, was set to 0.05. Let HN1,N2 indicate the Hamming distance and IM(N1,N2) indicate the Ipsen-Mikhailov distance between two undirected networks N1,N2, then the HIM metric of similarity HIM(N1,N2) is described by the following equation.
DARLIN data analysis and validation
Processed h5ad of the DARLIN dataset was downloaded from Zendo (tissue_adata_refined_20221106_joint.h5ad). The counts matrix in the anndata was log-normalized using scanpy. cellHarmony (community alignment as described above) was used to project our CITE-seq compendium cell atlas annotations on the cells from DARLIN dataset. The log-normalized gene expression matrix was provided as the input to cellHarmony.
Clone-by-cell count matrix was extracted from the above-mentioned h5ad (adata.obsm[‘X_clone’]). Using this matrix, log2-normalized clonal coupling scores were calculated using the approach described above used CellTag lineage barcoding data. Lineage relationships that had log-clonal coupling score greater than 0 and were evidenced by at least 3 clones were considered for visualization purposes (Fig. 7c, Extended Data 9d). Further, lineage relationships with log-clonal coupling score greater than 0, greater than 0.3, and greater than or equal to 0.3 were assigned as level 1, level 2, and level 2 (relaxed), respectively in Supplemental Table 7.
CITE-seq and TEA-seq Integration
To transfer labels from the titrated CITE-seq and matching TEA-seq captures, the two modalities were integrated using the software harmonypy, separately for each capture (HSC-MPP and MultiLin). Principal components (nPCs=30) were calculated using the sklearn.decomposition.PCA function after concatenating the CPTT-normalized gene expression matrices from CITE-seq and TEA-seq captures. Principal component 1 was removed as it was correlated strongly with CITE-seq vs. TEA-seq batches, and remaining PCs (PC2-PC30) were used as input to harmonypy to generate an embedding with minimal platform batch effects. A nearest neighbor classification was then carried out from CITE-seq defined clusters to TEA-seq cells (KNN=10) using the sklearn.neighbors.KNeighborsClassifier class, identifying of 57 corresponding clusters. UMAP embeddings were derived using a custom python script (Synapse: syn53237568).
Cluster Specific TEA-seq ATAC Processing
Using the cluster assignments from the CITE-seq to TEA-seq classification, we split the 10X Cell Ranger chromatin accessibility position-sorted BAM file into separate files for each cluster. Peaks were called using the MACS2 peak calling algorithm with the following “--nomodel --shift 37 --ext 73 -p 0.05 -B --SPMR --call-summits”. To merge peaks called between cluster splits, summits were ranked by p-value from MACS2 and overlapping summits with higher p-values were removed within a window of +/−500 bp, yielding approximately 800,000 peaks. To focus on dynamic peaks of interest, the log2-CPM normalized read count value for each peak was compared to the corresponding log2-CPTT CITE-seq gene expression values using a Pearson correlation test across the 57 clusters. These correlations were restricted to peaks with a given gene’s previously defined TAD identified in murine HSPCs^59^ (GSE119347, “BMHSC_TADs.bed.gz” translated from mm9 to mm10 using UCSC lift-over tool). Only the top 5,000 genes ranked by variance were considered. Any peak with a p-value < 0.001 from this test was nominated as a dynamic peak used for downstream analyses, yielding approximately 100,000 dynamic peaks.
Base Resolution Contribution Scores for Tn5 Insertion and CWM Pattern Scoring
To estimate the sequence specific importance for Tn5 insertion from chromatin accessibility captured by TEA-seq, bias corrected ChromBPNet models were generated on each pseudobulk ATAC bam split file for clusters with enough reads to complete model training (32 of the 57 mapped clusters). Bias models were created using the total possorted BAM file from one replicate of the MultiLin TEA-seq capture. The Tn5 counts prediction head model was then used to create contribution scores at single base resolution across the dynamic peak set. These contribution scores were scanned for seqlets of frequently used base pair patterns using the TF-MoDISco-lite program, yielding approximately 30–50 CWM patterns for each cluster (approximately 1,000 patterns total). To identify all seqlets across all dynamic peaks, as opposed to the sampled set used by TF-MoDISco for clustering, the base frequency of each CWM pattern was used as input to scan for loci with similar base pair patterns using the gimmemotifs python package. Matching loci were subsequently scored in their ability to match the CWM by taking the dot product of the base resolution contribution score at matching nucleotides of that loci to the corresponding positions in the CWM. Any loci that matched the CWM and had a dot product score greater than the 5^th^ percentile of those seqlets of the CWM previously called by TF-MoDISco-lite were included in the final set. Thus, this created a set of approximately 20 million seqlets, each matching a CWM pattern identified by TF-MoDISco-lite from the 32 cluster-specific ChromBPNet models. Further details and code are provided in our GitHub repository.
Constructing Putative Gene Regulatory Networks
To infer CWM patterns that likely contribute to gene expression, pairwise Pearson correlation tests were performed between all seqlets within a given gene’s TAD (as previously defined in dynamic peak selection under “Cluster Specific TEA-seq ATAC Processing”) to that gene across the 32 pseudobulk clusters for which ChromBPNet models could be trained. Only the top 5,000 genes ranked by variance (32 pseudobulks) were considered. Any seqlet that correlated to a gene within a TAD above a threshold of 0.4 was maintained as a candidate connection. For visualization, we further restricted to the top 20 correlated seqlets to each gene, yielding approximately 80,000 seqlets to 5,000 genes. Pearson correlation values were recalculated across both these sets to generate a pair-wise correlation matrix between genes and seqlets. Each seqlet was annotated using its matching CWM and the CIS-BP2 motif database to define candidate transcription factors and their families The genes were ordered using hierarchical clustering of seqlet correlation values using the SciPy Python package. Seqlets were then grouped using MarkerFinder with gene clusters as groups. Base resolution contribution score values were then visualized using the UCSC genome browser. For heatmap visualization, we plotted the z-score log2-CPTT normalized gene expression values with the same hierarchically clustered order of genes from the seqlet to gene expression correlation heatmap to show their expression across the 32 clusters. We replicated this procedure for the seqlet contribution score values for the MarkerFinder ordered 80,000 seqlets. Next, the contribution score values for the seqlets across the clusters were aggregated using the CIS-BP2-defined transcription factor families to visualize the average contribution scores for seqlets of a given transcription factor family. This procedure was repeated at the level of CWM patterns defined by TF-MoDISco-lite. To derive the shown GRN model, we restricted visualization to the most informative lineage markers for all 32 clusters, by intersecting the 5,000 top variably expressed genes with the top 25 MarkerFinder markers from all 88 CITE-seq clusters, yielding ~500 putative target genes. The seqlet correlation matrix was rederived with this 500 gene by 80,000 seqlet set and replotted as a pair-wise correlation heatmap and the Cytoscape network produced from all CWM associated transcription factors and top marker genes (object deposited in Synapse syn60529836).
In Silico Perturbation of Select Transcription Factors
To test the cluster specific GRNs built using ChromBPNet defined seqlets as well as transcription factor and target gene expression, we input these connections as the base GRN provided to CellOracle. The Oracle object was constructed using all expressed genes and an auto-selected k-value of 80 with all other parameters set to defaults. 100 cells were randomly sampled from each of the 32 clusters in the HSC-MPP and MultiLin gate populations for pruning and in silico perturbation steps. Cell transition shifts were calculated by simulating KO of select transcription factors (Gata1, Gata2, Irf8, Spi1, Cebpa, and Cebpe) by setting the normalized expression value to 0 in all cells and the n_propagation parameter set to 5. Following in silico perturbation, the shifted transcriptomes were tested against the input transcriptomes using an empirical Bayes adjusted linear model (limma^122^) to assess significance of change in target gene expression.
Validation of GRN Seqlet Loci
To quantify the overlap between GRN seqlet positions and previously defined enhancer regions (cCREs) we used bedtools intersect to count the overlap between seqlet positions and cCRE regions. We used GIGGLE index and search functions with default parameters to estimate the significance of overlap between the collection of seqlet positions associated to each transcription factor to the peaks that were previously defined in all murine ChIP-seq samples targeting transcription factors. The significance of self-to-self enrichment between a given transcription factors seqlet positions and the corresponding ChIP-seq sample regions was tested using a one-sided Mann Whitney U test of the self-to-self GIGGLE score against the self-to-others scores.
Index Sorting and Assessment of CRE Reporter Fate
Index sorting was carried out on a Sony MA900 or BD FACSymphony S6 sorter on mice with lineage-specific-cre systems that activate tdTomato upon lineage commitment. Single cells were sorted into 96-well plates under “in-vitro culturing conditions of HSPCs” (above) and cells were analyzed either under an Olympus fluorescent microscope after 5 days to count the number of cells produced and how many became tdTomato^+^, or alternatively using the Agilent BioTek Cytation C10.
N. brasiliensis Infection Model Time Course
750 infectious larvae of N. brasiliensis were inoculated by subcutaneous injection to 8-week old Balb/cJ male mice.^75,76^ Bone marrow was collected from mice at time points of 3, 5, 7, and 10-days post infection or from uninfected mice as control samples (day 0). Bone marrow was CD117-enriched as described above and processed by full spectrum flow cytometry with a Cytek Aurora cytometer and captured using the Chromium X with Chromium Next GEM Single Cell 3ʹ kit v3.1 chemistry (PN-1000268, 10X Genomics). Supervised assignment of cell labels to the multimodal scTriangulate annotations were performed using cellHarmony as described above. cellHarmony differential expression analyses were performed for each assigned cell state relative and time-point relative to day 0 (eBayes t-test p<0.05, FDR corrected and fold>1.2), with associated heatmap visualization in AltAnalyze.
Extended Data
