Reconstruction of complex single-cell trajectories using CellRouter

Abstract

A better understanding of the cell-fate transitions that occur in complex cellular ecosystems in normal development and disease could inform cell engineering efforts and lead to improved therapies. However, a major challenge is to simultaneously identify new cell states, and their transitions, to elucidate the gene expression dynamics governing cell-type diversification. Here, we present CellRouter, a multifaceted single-cell analysis platform that identifies complex cell-state transition trajectories by using flow networks to explore the subpopulation structure of multi-dimensional, single-cell omics data. We demonstrate its versatility by applying CellRouter to single-cell RNA sequencing data sets to reconstruct cell-state transition trajectories during hematopoietic stem and progenitor cell (HSPC) differentiation to the erythroid, myeloid and lymphoid lineages, as well as during re-specification of cell identity by cellular reprogramming of monocytes and B-cells to HSPCs. CellRouter opens previously undescribed paths for in-depth characterization of complex cellular ecosystems and establishment of enhanced cell engineering approaches.

Fig. 1

Overview of CellRouter. Step (1) Starting from single cells representing multiple cell states, including stable and in transition, a gene regulatory network was built to identify putative gene regulatory relationships. Step (2) Subpopulations were identified by a combination of learning the structure of the data and community detection algorithms. Step (3) High-resolution multi-state trajectories are identified from any subpopulation to any other, including intermediate and mature cell types. Step (4) Multi-state trajectories illustrate the divergence of single-cell transcriptomes along each trajectory progression. Identification of genes implicated in the dynamic biological process under study, such as differentiation, and identification of regulators driving or mediating cell-fate transitions at the gene and network level

Fig. 2

Gene expression dynamics during neutrophil differentiation. a t-SNE plot using genes reliably expressed as identified in the original study. b Predicted transcriptional regulators during cell-fate transitions from hematopoietic stem cells to neutrophil progenitors (subpopulation 17) and mature neutrophils (subpopulation 18). c Transcriptional dynamics from hematopoietic stem cells (subpopulation 20) to mature neutrophils (subpopulation 18). d Validation of the developmental timing of Cebpe and Mxd1 by bulk RNA-seq along a time-course of neutrophil differentiation. e t-SNE map colored by expression of Mxd1. f t-SNE maps colored by expression of Mxd1 predicted target genes during neutrophil differentiation. g Kinetic trends of top regulators of transitions from HSCs to subpopulation 18 in b, top panel, and validation of these patterns by bulk RNA-seq during a time-course of neutrophil differentiation, bottom panel. h Genes known to be important in progenitor and mature neutrophils, top panel, and the validation of their expression dynamics using the time-course of neutrophil differentiation, bottom panel. Accession codes: GSE76983 and GSE84874

Fig. 3

Gene expression dynamics during erythroid differentiation. a Predicted transcriptional regulators of transitions from HSCs to erythroblasts (subpopulations 20 to 10). b Transcriptional dynamics of top regulators of transitions from HSCs (subpopulation 20) to the erythroblast subpopulation 10. c Predicted transcriptional regulators of transitions from HSCs to the erythroblast progenitor subpopulation 19. d Staining and representative fluorescence-activated cell sorting (FACS) profiles of human CD34⁺ peripheral blood cells induced to differentiate to erythroid cells (scale bar: 100 μm). e Transcriptional dynamics along the CellRouter trajectory of genes involved in mouse and human erythroid differentiation from HSCs to erythroblasts (subpopulations 20 to 10, top panel) and the qPCR validation of these patterns during human erythroid differentiation, bottom panel. f Transition-specific gene expression, as shown by the presence of an expression pattern in the top panel and lack of any patterns in the bottom panel. g Gene ontology enrichment analysis of genes positively regulated during HSC differentiation (subpopulation 20) to erythroid progenitors and more mature erythroblasts (subpopulations 10 and 19) as well as neutrophil progenitors and mature neutrophils (subpopulations 17 and 18). Expression of each curve was normalized between 0 and 1 to highlight their dynamic patterns. Accession code: GSE76983

Fig. 4

Multi-lineage hematopoietic stem cell differentiation. a Three-dimensional diffusion components analysis with subpopulation structure identified by CellRouter. b Identification of lineages based on expression of known marker genes with Procr expressed in HSCs and Gata1, Ccl3, and Ctsg expressed in the erythroid, lymphoid, and granulocyte-macrophage lineages, respectively. c Gene ontology analysis using genes dynamically upregulated in selected differentiation trajectories from HSCs to the erythrocytes (subpopulation 8), granulocyte/macrophage progenitors (subpopulation 11), and lymphoid progenitors (subpopulation 4). d Predicted regulators of transitions to the erythroid lineage (subpopulation 8). e Predicted transcriptional regulators of lymphoid development from HSCs (left panel) and their pseudo-temporal regulation (right panel). f Clustering of gene expression trends during lymphoid development from HSCs into five expression patterns, with a summary of the patterns in each cluster shown in the bottom right panel. g Representative genes in transcriptional clusters from f and the Jaccard similarity matrix of the lymphoid trajectory demonstrating the gradual divergence of single-cell transcriptomes during differentiation (bottom right panel). Accession code: GSE81682

Fig. 5

Application of CellRouter to gain insight for cell reprogramming. a STEMNET dimensionality reduction with subpopulations identified by CellRouter and cell types annotated based on the original study. b Positive and negative controls for lineage-specific differentiation trajectories. c Left panel: predicted candidates to reprogram B-cells (subpopulation 10) to HSPCs (subpopulation 12) based on overexpression (genes in red) or knockdown (genes in blue) of indicated transcriptional regulators. Middle panel: Derivative of the gene expression dynamics for each regulator. Right panel: Kinetic profile of each regulator along the reprogramming trajectory from B-cells to HSPCs. d Left panel: predicted candidates to reprogram mono/DCs (subpopulation 17) to HSPCs (subpopulation 12) based on overexpression (genes in red) or knockdown (genes in blue) of indicated transcriptional regulators. Middle panel: Derivative of the gene expression dynamics for each regulator. Right panel: Kinetic profile of each regulator along the reprogramming trajectory from mono/DCs to HSPCs. e Regulator-centered subnetworks. * indicates genes previously used to convert or reprogram a starting population to blood progenitors or HSPCs. Accession code: GSE75478

Fig. 6

Benchmarking CellRouter to other algorithms. a Simulated curves representing gene expression dynamics of a selected gene along differentiation trajectories demonstrating the relationship between correlation and autocorrelation with the trajectory. b t-SNE map of myeloid progenitor cells annotated by subpopulations identified by CellRouter. Inset shows broad cell type annotations for each population. c Positive and negative controls of erythrocyte and GMP differentiation. d Number of genes up- or downregulated significantly correlated with the trajectories identified by each trajectory detection algorithm. e Left panel: selected markers of granulocyte/monocyte progenitors and erythrocytes. Size is proportional to the percent of cells expressing the markers in each subpopulation. Color is proportional to mean expression of the markers in each subpopulation. Right panel: lag-1 autocorrelation of marker genes on the right panel along the GMP and erythrocyte differentiation trajectories identified by each trajectory detection algorithm. Accession code: GSE72857