Herein, we briefly describe the main steps involved in the calculation of single-cell transcriptional uncertainty using CALISTA ( Papili Gao et al., 2020).

Pre-processing. Given an N × G single-cell expression matrix M , where N denotes the number of cells and G the number of genes, the pre-processing in CALISTA involves two steps: a normalization of the expression data m n , g – i.e. the number of transcripts of gene g in the n-th cell, and a selection of the most variable genes ( Papili Gao et al., 2020).

Cell clustering. CALISTA clustering follows a two-step procedure. The first step involves a greedy optimization strategy to find cell clustering that maximizes the total cell likelihood, i.e. the sum of the likelihood value for all cells. The single-cell likelihood value is computed as the joint probability of the cell’s gene expression data, which is set equal to the product of the probabilities of the mRNA counts for the selected genes based on the mRNA distribution from the two-state stochastic gene transcription model. To avoid issues with numerical overflow, we use the logarithm of the cell likelihood. By performing the greedy optimization multiple times, a consensus matrix containing the number of times two cells in the dataset are put in the same cluster, is generated. In the second and final step, CALISTA generates the cell cluster assignments by using k-medoids clustering based on the consensus matrix. The final outcome of CALISTA’s clustering is the assignment of cells into K clusters and the optimal model parameters for the two-state gene transcription model: θ g k = θ on θ off θ t g k . for each gene g in cluster k ( Papili Gao et al., 2020). In this case, θ on is the normalized rate of the promoter activation, θ off is the normalized rate of the promoter inactivation, and θ t is the normalized rate of mRNA production when the promoter is active. These parameters are normalized by the rate constant of mRNA degradation θ d , so that θ d = 1 .

Lineage progression inference. In CALISTA, cell lineage progression is inferred based on cluster distances – a measure of dissimilarity between two clusters. The cluster distance of two cell clusters is defined as the average decrease in the cell likelihood value if the cells from these two clusters are grouped as one cluster, as opposed to the original clustering. The lineage progression graph is built by adding transition edges between pairs of clusters in increasing magnitude of cluster distance until all clusters are connected to at least one other cluster, or based on user-specified criteria.

Single-cell transcriptional uncertainty. The last step in our analysis is to compute the final single-cell likelihood. Briefly, for each cell, we consider all edges in the lineage progression graph that are adjacent to the cell’s respective cluster, i.e. edges that eminate from or pointing to the cluster to which the cell belongs. The likelihood of a cell along an edge is evaluated by interpolating the likelihood values of the cell’s gene expression using the mRNA distributions from the two adjacent clusters. Each cell is then assigned to the edge along which its interpolated likelihood value is maximum, and the final cell likelihood is set to this maximum value. As mentioned above, the single-cell transcriptional uncertainty is evaluated as the negative logarithm of the cell likelihood value (NLL).

Pseudotimes calculation. We can evaluate the pseudotimes for the cells according to the following procedure. First, a pseudotime is given to each cluster with a value between 0 (initial cell state) and 1 (final cell fate). Subsequently, we determine the linear fractional position of each cell along its respective edge at which its interpolated likelihood value is maximum (see Single-cell transcriptional uncertainty). The pseudotime of a cell is computed by a linear interpolation of the pseudotimes of the two clusters adjacent to its assigned edge according to the cell’s linear fractional position on this edge.

Epigenetic landscape reconstruction. To visualize the 3D transcriptional uncertainty landscape, we apply dimensional reduction techniques such as principal component analysis (PCA) or t-SNE on the z-scored expression data, to project the gene expression of each individual cell on two dimensional axis, which gives the x-y axis of the landscape plot. For the z axis, we plot the NLL values. The transcriptional uncertainty landscape surface is reconstructed by estimating local approximation of individual cell 3D coordinates on a regular 30×30 grid by using a publicly available Matlab (R2020a) surface fitting package called gridfit.

Bargaje et al. scRT-qPCR dataset. The dataset includes the expression profiles of 96 genes from 1896 single cells at eight different time points (day 0, 1, 1.5, 2, 2.5, 3, 4, 5) during the differentiation of human pluripotent stem cells (iPSCs) into either mesodermal (M) or endodermal (En) fate ( Bargaje et al., 2017). By employing CALISTA, we obtained five cell clusters and detected a bifurcation event, which gives rise to the two final cell fates. After lineage inference, we pseudotemporally ordered cells along the inferred differentiation paths (for more details, see ( Papili Gao et al., 2020)).

Treutlein et al. scRNA-sequencing dataset. The dataset includes the gene expression profiles of 405 cells during reprogramming of mouse embryonic fibroblast (MEF) into a desired induced neural (iN) and an alternative myogenic (M) cell fate ( Treutlein et al., 2016). We pre-processed the data using CALISTA to select the 40 most variable genes (10% of the number of cells) for the transcriptional uncertainty analysis. CALISTA identified four different subpopulations and successfully recovered the bifurcation event (for more details, see ( Papili Gao et al., 2020)).

Richard et al. scRT-qPCR dataset. The dataset contains the expression profile of 91 genes measured from 389 cells at six distinct time points (0, 8, 24, 33, 48, 72 h) during the differentiation of primary chicken erythrocytic progenitor cells (T2EC) ( Richard et al., 2016). Following the CALISTA pre-processing step, we removed cells in which less than 75% of the genes are expressed. Then, we selected the subset of genes with at least one non-zero expression values. A total of 354 cells and 88 genes were considered in the transcriptional uncertainty analysis. Based on eigengap heuristics ( Papili Gao et al., 2020; von Luxburg, 2007), we grouped cells into six optimal clusters and ordered cells along the inferred linear trajectory (see the Extended data S8 Figure ( Gao et al., 2023)).

Stumpf et al. scRT-qPCR dataset. The dataset comprises the single-cell expression of 97 genes at seven time points (0, 24, 48, 72, 96, 120, 168 h) during neural differentiation of mouse embryonic stem cells (E14 cell line) ( Stumpf et al., 2017). In the data pre-processing, we excluded cells in which less than 70% of genes are expressed. Then, we selected genes with at least one non-zero expression values. A total of 276 cells and 93 genes were considered for for the transcriptional uncertainty analysis. Based on eigengap heuristics ( Papili Gao et al., 2020), we grouped cells into five optimal clusters and ordered cells along the inferred linear trajectory ( Extended data S9 Figure ( Gao et al., 2023)).

Moussy et al. scRT-qPCR dataset. The single-cell expression dataset includes normalized Ct values for 91 genes in 435 cells captured at five distinct time points (0, 24, 48, 72, 96 h) during human cord blood-derived CD34+ differentiation ( Moussy et al., 2017). We employed CALISTA to group cells into seven clusters, reconstruct the developmental trajectory and calculate pseudotimes ( Extended data, S10 Figure ( Gao et al., 2023)).

Guo et al. scRT-qPCR dataset. The dataset comprises the single-cell expression values of 48 genes from 387 individual cells isolated at four distinct developmental cell stages, from 8-cell stage mouse embryos to 64-blastocyst ( Guo et al., 2010). By applying CALISTA, we identified seven different subpopulations along the differentiation process, and the inferred lineage hierarchy pinpointed two bifurcations events at 32- and 64-cell stage ( Extended data, S11 Figure ( Gao et al., 2023)). The timing of the lineage bifurcations coincides with two well-known branching points: one at 32-cell stage when totipotent cells differentiate into trophectoderm (TE) and inner cell mass (ICM), and another at 64-cell stage when ICM cells differentiate into primitive endoderm (PE) and epiblast (E).

Nestorowa et al. scRNA-sequencing dataset. The dataset comprises single-cell gene expression of 1656 cells from mouse hematopoietic stem cell differentiation ( Nestorowa et al., 2016). We pre-processed the data by removing genes with non-zero values in less than 10% of the cells. Then, we selected 433 most variable genes, which is 10% of the number of genes after the previous pre-processing step, for the transcriptional uncertainty analysis ( Papili Gao et al., 2020). We set the optimal number of clusters based on the original study ( Nestorowa et al., 2016), which reported six different subpopulations and two bifurcation events: the first one producing common myeloid progenitor (CMP) from lymphoid-primed multipotent progenitors (LMPP), and the second one generating granulocyte–monocyte progenitors (GMP) from megakaryocyte-erythroid progenitors (MEP) ( Extended data, S12 Figure ( Gao et al., 2023)).

Moignard et al. scRT-qPCR dataset. The dataset contains the single-cell expression level of 18 transcription factors measured in a total of 597 mouse bone marrow cells during hematopoietic differentiation. By applying CALISTA, we successfully identified the five subpopulations and the two branching points detected in the original study ( Moignard et al., 2013): long-term hematopoietic stem cells (HSC) differentiating into megakaryocyte–erythroid progenitors (PreM) or lymphoid-primed multipotent progenitors (LMPP); LMPP cells differentiating into granulocyte–monocyte progenitors (GMP) and common lymphoid progenitors (CLP) (for details see ( Papili Gao et al., 2020)).

We defined gene transcriptional burst size and burst frequency using the two-state model parameters, as follows: Burst size = θ t θ off (1) Burst frequency = θ on (2)

The burst sizes and burst frequencies were evaluated using the cluster parameters θ g k = θ on θ off θ t g k obtained from single-cell clustering analysis of CALISTA. Meanwhile, the average gene-wise NLL values for each single-cell cluster was computed as: NLL g k = ∑ n = 1 N k NLL g n N k (3)

where NLL g n is the negative log-likelihood of cell n based only on the expression of gene g, and N _k is the total number of cells in cluster k.

Cells and genes were first filtered based on the pre-processing strategy in the original publication by La Manno and colleagues ( La Manno et al., 2018), which resulted in a total of 1720 cells and 1448 genes from human glutamatergic neurogenesis, and a total of 18140 cells and 2141 genes from the mouse hippocampus dataset. We further reduced the number of genes to only the top 500 highly variable genes for the transcriptional uncertainty analysis. The cell cluster assignments generated by Velocyto – the algorithm for computing RNA velocity from the original publication ( La Manno et al., 2018) – were considered, instead of using CALISTA. Based on the clustering, we employed CALISTA to generate the lineage progression and cell pseudotimes ( Extended data, S13 Figure ( Gao et al., 2023)). The RNA velocity and transcriptional uncertainty values for the top 500 genes were calculated by employing Velocyto and CALISTA, respectively. The cell-wise RNA velocity was set to the Euclidean norm of the vector of RNA velocities for each cell, while the cell-wise NLLs was NLL n k = ∑ g = 1 500 NLL g n 500 (4)

In this work, we used CALISTA ( Papili Gao et al., 2020), a likelihood-based bioinformatics toolbox designed for an end-to-end analysis of single-cell gene expression data, to evaluate the transcriptional uncertainty of each individual cell based on its gene expression data (see the Extended data, Supplementary Notes S1 ( Gao et al., 2023)). CALISTA uses the two-state model of stochastic gene transcription bursts to characterize the steady state distribution of mRNA counts in individual cells ( Peccoud & Ycart, 1995). In the model, a gene promoter stochastically switches between the ON and OFF state, and only in the ON state can gene transcription occur. The distribution of mRNA depends on four model parameters: θ _on (the rate of promoter activation), θ _off (the rate of promoter inactivation), θ _t (the rate of mRNA production when the promoter is in the ON state), and θ _d (the rate constant of mRNA degradation) ( Herbach et al., 2017; Kim & Marioni, 2013) (see Figure 1a). For example, when θ _off >> θ _on and θ _off >> θ _d , keeping θ _t / θ _off fixed, mRNA are produced through bursts of short but intense transcription, which is a typical case observed for gene transcriptions in single cells ( Munsky et al., 2012) ( Suter et al., 2011). As the mRNA distribution is linked to mechanistically interpretable parameters, CALISTA is able to give insights into the possible mechanism driving the cell heterogeneity dynamics during cell differentiation.

CALISTA employs a maximum likelihood approach and assigns a likelihood value to each cell based on its gene expression based on the mRNA distribution governed by the two-state model of stochastic gene transcription. The single-cell transcriptional uncertainty is evaluated as the negative logarithm of the likelihood (NLL) value for a cell. The single-cell likelihood value reflects the joint probability of its gene expression repertoire. A cell with a low likelihood value may indicate that the gene expression of the cell is different from its neighboring cells, i.e. the cell is an outlier. But, more interestingly, a low likelihood value may also correspond to a cell state of high uncertainty in the gene expression. The group of cells in such a high uncertainty state have gene expressions that are dissimilar to each other, and thus, the gene expression distribution will have a high entropy. By visualizing the single-cell transcriptional uncertainty over the two-dimensional projection of the single-cell transcriptomics data—for example, using the first two principal components from PCA—we constructed a transcriptional uncertainty landscape in the form of a surface plot of the NLL value. In this way, we studied the landscape of transcriptional uncertainty during cell differentiation at single-cell resolution. On such single-cell transcriptional uncertainty surface, an aberrant cell can be easily distinguished from a cell of high uncertainty state, since such an aberrant cell will appear isolated from its neighboring cells with a high NLL.

In the following, we demonstrated an application of our procedure described above to a single-cell transcriptional dataset from cardiomyocytes differentiation from human induced pluripotent stem cells (iPSCs) ( Bargaje et al., 2017). The single-cell clustering of CALISTA returned five clusters ( Papili Gao et al., 2020) and identified one bifurcation event in the lineage progression, which led to two cell lineages ( Bargaje et al., 2017), in good agreement with the number of cell types reported in the original study. The estimated uncertainty landscape shows cells exiting the initial epiblast state that is characterized by a valley in the landscape, passing through a hill of high transcriptional uncertainty corresponding to primitive streak (PS)-like progenitor state, before ending up at one of the low transcriptional uncertainty terminal states corresponding to either mesodermal (desired) or endodermal (undesired) fate (see the Extended data, S1 Figure ( Gao et al., 2023)). As depicted in Figure 1b, the intermediate cell cluster (cluster 2) comprising PS-like cells have higher cell uncertainty (lower single-cell likelihood) than the other clusters. Figure 1c and d give the moving-averaged uncertainty values for pseudotemporally ordered cells using a moving window of 10% of the total cells for both endodermal and mesodermal paths, respectively. The moving-averaged transcriptional uncertainty for the two differentiation paths follows a rise-then-fall trajectory where the peak of uncertainty coincides with the lineage bifurcation event.

We explored whether the rise-then-fall in uncertainty is an artefact from using the two-state model to evaluate the cell likelihood values. To this end, we implemented a modified version of the algorithm for ordering cells by calculating the cell likelihood values using the empirical (observed) distribution, instead of the analytical distribution from the two-state model. As shown in the Extended data S2 Figure ( Gao et al., 2023), the transcriptional uncertainty landscape from the modified implementation shows a strong resemblance to the original one. We also investigated whether the number of clusters may affect the landscape, in which using too few of the clusters may artificially inflate the uncertainty due to the mixing of cells from different states. We reran CALISTA by using a higher number of clusters (set to nine based on the eigengap heuristic ( von Luxburg, 2007)). The hill in the uncertainty landscape is again seen around the bifurcation event upon using a higher number of cell clusters ( Extended data, S3 Figure ( Gao et al., 2023)). Finally, we used a different algorithm to cluster cells, specifically using a Laplacian-based clustering algorithm called single-cell interpretation via multikernel learning (SIMLR) ( Wang et al., 2017), to test whether the shape of the transcriptional uncertainty landscape changes with the clustering algorithm. The single-cell clusters can be interpreted as the transitional states that the differentiating cells go through. Starting with the result of SIMLR cell clustering, we then generated the lineage progression and estimated the cell likelihood values using CALISTA. The transcriptional uncertainty landscape from SIMLR cell clustering has the same shape as that in Extended data S1 Figure ( Gao et al., 2023), demonstrating that the transcriptional uncertainty landscape observed above is not dependent on using CALISTA for cell clustering ( Extended data, S4 Figure ( Gao et al., 2023)).

To further elucidate the role of specific genes in shaping the transcriptional uncertainty landscape, we looked at the transcriptional uncertainty associated with individual genes. Figure 1e depicts the NLL distribution of each gene for the five single-cell clusters. As expected, cells in cluster 2 have generally higher NLL than those in the other clusters. Figure 1e clearly illustrates that within cluster 2, some genes show higher NLL values than the others ( Extended data, S5 Figure ( Gao et al., 2023)). To identify the important genes related to transcriptional uncertainty, we identified genes with NLL values exceeding a threshold δ for at least 30% of the cells in each cluster, where δ is set to 3 standard deviation above the overall mean NLL for all cells and genes in the dataset (see Methods equation (2)). None of the genes in clusters 1, 4 and 5 have a NLL above the threshold. Meanwhile, 16 and eight genes in clusters 2 and 3, respectively, pass the above criterion for high uncertainty with four common genes between the two gene sets ( Extended data, S1 Table ( Gao et al., 2023)). Genes with high transcriptional uncertainty in cluster 2 may have functional roles in cell fate determination. The gene set of cluster 2 includes known genes upregulated only in the PS-like state (e.g. EOMES, GSC, MESP1 and MIXL1), as well as markers of mesodermal and endodermal cells (e.g. BMP4, HAND1, and SOX17) ( Bargaje et al., 2017; Papili Gao et al., 2020) ( Extended data, S6 Figure ( Gao et al., 2023)). Meanwhile, the main contributors to cell uncertainty in cluster 3 (e.g. BMP4 and MYL4 ( Bargaje et al., 2017; Papili Gao et al., 2020)) are known transition genes between PS-like cells and the final mesoderm fate ( Extended data, S7 Figure ( Gao et al., 2023)). Figure 1f depicts the protein-protein interaction (PPI) network related to the gene set of cluster 2 using STRING (minimum required interaction score of 0.4) ( Szklarczyk et al., 2015), indicating that these genes form a strongly interconnected hub of known transcription factors and molecules involved in the signal transduction of embryonic development ( Extended data, Table S1 ( Gao et al., 2023)).

We further applied the procedure above to seven additional single-cell transcriptomic datasets that were generated using scRT-qPCR ( Guo et al., 2010; Moignard et al., 2013; Moussy et al., 2017; Richard et al., 2016; Stumpf et al., 2017) and scRNA-sequencing ( Nestorowa et al., 2016; Treutlein et al., 2016), to assess the universality of the rise-then-fall feature of single-cell transcriptional uncertainty landscape during cell differentiation. The first of these datasets came from 405 cells during mouse embryonic fibroblast (MEF) reprogramming into induced neural (iN) and myogenic (M) cells ( Treutlein et al., 2016). Like the iPSC differentiation above, the lineage progression has a single bifurcation point. As depicted in Figure 2a, the single-cell transcriptional uncertainty increases from the initial MEF state and reaches a peak around the bifurcation before decreasing toward two end-point cell fates. The rise-then-fall of transcriptional uncertainty in the MEF reprogramming is in good agreement with what we observed in the iPSCs differentiation above. Higher entropy of gene expression distribution in a cell population has also been reported in the reprogramming of iPSCs ( Buganim et al., 2012).

Next, we analyzed datasets from cell differentiation processes without a lineage bifurcation and with multiple lineage bifurcations. Three scRT-qPCR datasets came from differentiation systems without bifurcation, including the Richard et al. study on chicken erythrocytic differentiation of T2EC cells ( Richard et al., 2016), the Stumpf et al. study on differentiation of mouse embryonic stem cells (ESC) to neural progenitor cells (NPC) ( Stumpf et al., 2017), and the Moussy et al. study during CD34+ cell differentiation ( Moussy et al., 2017). The single-cell clustering and lineage progression by CALISTA produced the expected cell differentiation trajectory (see Extended data, S8 to S10 Figures ( Gao et al., 2023)). The single-cell transcriptional uncertainty landscapes of these three differentiation systems, as shown in Figure 2b, exhibit a rise-then-fall profile, creating a hill that the cells traverse through in the differentiation process. A transitory increase in single-cell gene expression uncertainty was reported either directly or indirectly in the original publications. In Richard et al. (2016) and Stumpf et al. (2017), the authors adopted the Shannon entropy to quantify cell-to-cell variability (uncertainty), while Moussy et al. (2017) reported an unstable transition state with ‘hesitant cells’ flipping their morphology between polarized and round shapes before committing to the common myeloid progenitors-like fate. Morphological uncertainty therefore corresponded to a higher transcriptional uncertainty. Note that the Moussy et al. study looked at only the initial phase of the (hematopoietic) cell differentiation, and thus, it is likely that the differentiation process had not completed for the cells in the dataset.

The next set of single-cell gene expression data came from differentiation systems with multi-branching lineage, including the Guo et al. study during mouse embryo development from zygote to blastocyst ( Guo et al., 2010), Nestorowa et al. (2016) and Moignard et al. (2013) studies on hematopoietic stem cell differentiation. Figure 2c shows the single-cell transcriptional landscape for each of the datasets. For the Guo et al. study, we identified seven cell clusters and identified two bifurcations in the lineage. Here, we observed two hills in the transcriptional uncertainty landscape, each coinciding with a bifurcation event in the lineage progression – one at 32-cell stage (cluster 2 to cluster 3 and 4) and another at 64-cell stage (cluster 4 to cluster 6 and 7) (see the Extended data, S11 Figure ( Gao et al., 2023)). For the Nestorowa et al. (2016) ( Extended data, S12 Figure ( Gao et al., 2023)) and Moignard et al. (2013) (see Methods and ( Papili Gao et al., 2020)) datasets, we again observed peaks in the transcriptional uncertainty landscape that colocalize with the bifurcation points in the lineage progression.

The use of the two-state mechanistic gene transcriptional model within CALISTA enabled us to probe into a mechanistic explanation for the observed shape of the transcriptional uncertainty landscape. Table 1 show the pairwise Pearson correlations between the cell-averaged NLL of each cluster with two biologically interpretable model parameters, namely transcriptional burst size (number of transcripts generated in each burst) and burst frequency (occurrence of burst per unit time) ( Nicolas et al., 2017) (see Methods). The Pearson correlations indicate that the single-cell gene expression uncertainty increases with higher burst size and burst frequency ( p-value ≤ 0.01). Higher transcriptional burst size and frequency are associated with a lower θ _off – a lower rate of promoter turning off – and a greater θ _on – higher rate of promoter turning on. One possible explanation for such a change in model parameters is a higher chromatin accessibility during the transition period of cell differentiation. This finding is consistent with the view that stem cells increase its gene expression uncertainty or stochasticity by adopting a more open chromatin state to enable the exploration of the gene expression space ( Antolović et al., 2017; Fritzsch et al., 2018; Nicolas et al., 2017; Zhang & Zhou, 2018).

In a recent paper ( La Manno et al., 2018), La Manno and colleagues introduced the concept of RNA velocity, which involves computing the rate of change of mRNA from the ratio of unspliced to spliced mRNA. A positive RNA velocity indicates an induction of gene expression, while a negative RNA velocity indicates a repression of gene expression. La Manno et al. demonstrated that RNA velocities are able to predict the trajectory of cells undergoing a dynamical transition, such as in circadian rhythms or cell differentiation. In the following, we explored the relationship between RNA velocities and single-cell transcriptional uncertainty.

We evaluated the single-cell transcriptional uncertainty and RNA velocity for two single-cell gene expression datasets that were previously analyzed in La Manno et al. (2018). The first dataset came from human glutamatergic neurogenesis which has a linear (non-bifurcating) lineage progression. Figure 3 (top row) depicts the cell clustering, single-cell transcriptional uncertainty, and RNA velocities (see also Extended data, S13 Figure ( Gao et al., 2023)). The single-cell transcriptional uncertainty landscape again has the rise-then-fall shape, as in the other cell differentiation systems discussed above. Interestingly, the same rise-then-fall profile is also seen in the RNA velocity. As illustrated in Figure 3, the increase and decrease of the RNA velocity preceed the transcriptional uncertainty, and the peak of RNA velocity occurs prior to those of the transcriptional uncertainty (see the Extended data S1 File for animated illustration ( Gao et al., 2023)). Furthermore, a gene-wise cross-correlation analysis confirms a positive correlation between RNA velocity and single-cell transcriptional uncertainty with a delay for individual genes (see the Extended data, Figure S14 ( Gao et al., 2023)).

We also compared RNA velocity and single-cell transcriptional uncertainty for another dataset from mouse hippocampal neurogenesis with a multi-branching lineage ( La Manno et al., 2018). Figure 3 (bottom row) shows that like in the neurogenesis dataset earlier, the RNA velocity increases and then decreases during cell differentiation, and the change in the RNA precede that of the transcriptional uncertainty (see the Extended data, S2 File for animated illustration ( Gao et al., 2023)). Also, the RNA velocity peaks take place before the transcriptional uncertainty peaks. The rise-then-fall dynamic of the RNA velocity seen in the two datasets above is consistent with the view that cells engage in exploratory stochastic dynamics as they leave the progenitor state, and disengage this explorative mode as they reach toward the final cell state.