Discovering sparse transcription factor codes for cell states and state transitions during development

Computational analysis of gene expression to determine both the sequence of lineage choices made by multipotent cells and to identify the genes influencing these decisions is challenging. Here we discover a pattern in the expression levels of a sparse subset of genes among cell types in B- and T-cell developmental lineages that correlates with developmental topologies. We develop a statistical framework using this pattern to simultaneously infer lineage transitions and the genes that determine these relationships. We use this technique to reconstruct the early hematopoietic and intestinal developmental trees. We extend this framework to analyze single-cell RNA-seq data from early human cortical development, inferring a neocortical-hindbrain split in early progenitor cells and the key genes that could control this lineage decision. Our work allows us to simultaneously infer both the identity and lineage of cell types as well as a small set of key genes whose expression patterns reflect these relationships. DOI: http://dx.doi.org/10.7554/eLife.20488.001


Introduction
During development, pluripotent cells make a series of lineage decisions to give rise to the different cell types of the body. These lineage decisions are controlled by intra-cellular molecular networks that include transcription factors and signaling molecules. There are two fundamental challenges associated with understanding the differentiation of individual cells. The first is to identify lineage relationships: how cells and their progeny move from pluripotent through intermediate to terminally differentiated cell states. The second is to identify the key molecular drivers that allow cells to make fate decisions along their developmental trajectory.
Reconstructing cell lineages has traditionally involved prospectively tracking cells and their progeny using a variety of imaging or genetic tools (Buckingham and Meilhac, 2011;Frumkin et al., 2008;Orkin and Zon, 2008;Sulston et al., 1983). Recent progress in single-cell sequencing techniques (Grün et al., 2015;Jaitin et al., 2014;Macosko et al., 2015;Patel et al., 2014;Paul et al., 2015;Treutlein et al., 2014;Zeisel et al., 2015) allows for a complementary view of the transcriptional states of individual cells during the course of development, providing static snapshots of the dynamics of the underlying molecular network. But inferring lineage relationships and the dynamics of the underlying molecular networks has proved difficult using transcriptional data alone, in part because of the high dimensional nature of these data. Overcoming this challenge would be particularly useful for understanding the development of human organs such as the brain where traditional lineage tracing experiments are more difficult.
High dimensional data analysis techniques such as PCA, ICA, or t-SNE, are useful at reducing dimensionality. However, since the resulting axes represent linear combinations of a large number of features (for example, expression levels of each gene), interpreting the analysis or making experimental predictions is sometimes challenging. Meanwhile, traditional statistical methods such as linear multivariate regression have limited applicability for detecting patterns in high dimensional data (Advani and Ganguli, 2016;Donoho and Tanner, 2009). The challenges inherent in high-dimensional data analysis such as identifying discriminatory features are further exacerbated as the fraction of relevant features decreases (Donoho and Tanner, 2009). Computational techniques currently in use to cluster single cell data or to infer relationships among cells are built on these approaches, thereby assuming that all high-variance genes are equally relevant for pattern detection (Marco et al., 2014;Satija et al., 2015;Trapnell et al., 2014). In contrast, decades of work in developmental biology have revealed that combinations of a few transcription factors can be sufficient to experimentally perturb cell fate and developmental decisions (Gilbert, 2014;Graf and Enver, 2009;Takahashi and Yamanaka, 2006) suggesting that the expression patterns of a few genes may be most relevant for making computational inferences. Therefore, there is a need to detect patterns involving a small fraction of all genes. Unfortunately, except in the case of well-studied lineage decisions, we do not know the identity of this fraction.
In statistics, techniques relying on L1 regularization have been successful in contexts where the number of informative variables is known to be small but whose identities are unknown, both for regression problems (Baraniuk, 2007;Candès et al., 2006;Tibshirani, 1996;Wainwright, 2009) and for clustering (Witten and Tibshirani, 2010). Inspired by these successes in statistics, our aim here is to discover generalizable sparse patterns in gene expression data during development (if they exist), and to exploit these patterns to computationally infer the dynamics of cell state transitions from high-dimensional transcriptional data obtained during the course of development.
In this manuscript we analyze expression patterns among cell types with known lineage relationships in late hematopoiesis and discover a pattern in a sparse subset of genes that correlates with these relationships. We develop a Bayesian framework based on this gene expression pattern to simultaneously infer lineage transitions and the key genes that drive them. We apply this method to reconstruct the lineage tree among a different set of cell types in early hematopoietic development, and in this process identify many known drivers of early hematopoiesis, including Gata1, Cebpa and Ebf1. We further extend our method to analyze single-cell gene expression data, using genes exhibiting the discovered pattern to cluster cells from early brain development and to infer lineage relationships between these clusters. Our analysis reveals a split from early progenitors to putative neocortex and mid/hindbrain cell types, as evidenced by the mutually exclusive expression of region-specific genes such as FOXG1, LHX1, and POU3F2 (BRN2). This prediction was validated experimentally in a separate work (Yao et al., 2017). We finally discuss the advantages of using sparse patterns for making inferences and for modeling the underlying gene regulatory networks.

Discovering sparse patterns correlated with lineage transitions
In order to identify gene expression patterns that are robustly predictive of lineage relationships, we analyzed gene expression data from 41 cell types during B-and T-cell development that have an experimentally established developmental lineage ( Figure 1A, Heng et al., 2008). We searched for sparse patterns of gene expression amongst groups of three cell types from this collection; subsets of three are the minimal set in which measures of relative similarity can be used infer relative lineage relationships.
We identified 150 triplets of cell types with experimentally verified lineage relationships from Band T-cell development (Heng et al., 2008) (Figure 1-source data 1). Three such triplets are shown in Figure 1A. These triplets constituted both cell fate decisions (for example, cell type A gives rise to cell type B and C) and lineage progressions (cell type B gives rise to cell type A which then gives rise to cell type C). For each triplet, we noted which cell type was the progenitor or intermediate cell type ('root' cell type A) and which cell types were not ('leaf' cell types B and C). Note   (Heng et al., 2008). Three triplets -the minimal subset of the tree from which relative distances can be studied -are denoted, each including an intermediate root cell type (red) and terminal leaf cell types (blue). 150 triplets among all sets of three cell types within five steps on the lineage tree were extracted for pattern-detection. (B) For each triplet of cell types (left), each gene's expression level can have the clear minimum pattern in either the root or the leaves (right, top box) where the distribution of gene expression levels in one cell type is well separated from the other two (left, p<0.005 in a two sample t-test); clear maximum pattern (left, bottom) in the root or the leaves: the gene has a clear maximum in one of the cell types (right, p<0.005 in a two sample t-test). (C) A histogram of the number of genes showing the clear minimum pattern among the 150 triplets with known developmental topology. Triplets in which the root has the most genes showing the pattern shown in red; triplets in which one of the leaves has the most genes showing the pattern shown in blue. None of the triplets with more than 10 genes showing the pattern have the most genes with a clear minimum in the root (no red in any histogram bar except for the Figure 1 continued on next page that even in the case in which cell type B gives rise to cell type A which gives rise to cell type C, we will refer to A as the 'root' and B and C as 'leaves' for that particular triplet. We analyzed transcription factor gene expression data for these triplets from the Immunological Genome Consortium (Heng et al., 2008) since transcription factors are the ultimate drivers of cell fate decisions.
Surprisingly, we found that specific expression patterns involving only single genes correlated well with lineage relationships between three cell types. Genes that are not differentially expressed within a triplet of cell types convey no information about relationships between the cell types. Therefore, we select genes with expression variability among the three cell types. The expression pattern of such genes can belong to one of only two possible patterns: ( Figure 1B): the clear minimum pattern: the gene has a clear minimum level in one of the three cell types, with its distribution of gene expression levels being well-separated from the other two (p<0.005 in both two sample t-tests between the minimum cell type and the other two cell types; Triplets of cell types can be separated into categories based on how many genes exhibit the aforementioned patterns. 56% of the triplets contained more than 10 genes exhibiting the clear minimum pattern, and in 100% of these triplets the majority of genes with expression fitting this pattern reached their minimum expression in one of the leaves ( Figure 1C) and never in the root of the triplet. The frequency with which the pattern correctly indicated the lineage relationship increased with the number of genes within a triplet exhibiting the pattern, thus suggesting a confidence measure. The genes showing the clear minimum pattern fell into two distinct groups, corresponding to whether the minimum expression level was in one or the other of the two leaves ( The clear maximum pattern was a poorer indicator of lineage relationships (Figure 1-figure supplement 1C). 83% of the triplets had more than 10 genes exhibiting this pattern, but 10% of those showed the majority of genes with expression fitting this pattern reaching their maximum in the root while the others did so in the leaves. Crucially, the integrity of the relationship between the clear maximum pattern and lineage topology was not correlated with the number of genes exhibiting the pattern (Figure 1-figure supplement 1C). While the clear maximum pattern did not correlate with lineage relationships, genes exhibiting this pattern identify individual cell types, and therefore we will refer to them as marker genes.
There are many examples of genes known to be functionally important for lineage decisions whose expression patterns fit the clear minimum pattern. In the case of lateral inhibition commonly used during development, progenitor cells express genes together (for example, Notch and Delta) which are differentially expressed in the differentiated states (only Notch or only Delta) (Perrimon et al., 2012) reaching the minimum expression level in one of the leaves. The same pattern is also seen in multiple examples of lineage decisions often involving mutual inhibition, where key genes expressed in the progenitor are differentially regulated in the progeny (Graf and Enver, 2009;Qi et al., 2013;Thomson et al., 2011;Zhang et al., 1999). In each of these cases, key genes reach minimal expression levels in one of the leaves of the triplets.
The observation that genes exhibiting the clear minimum pattern are correlated with the lineage topology of a triplet of cell types further revealed that (i) only this fraction of transcription factors can be useful for inferring lineage relationships, and (ii) the identity of this fraction depends on which group of three cell types were analyzed. As the subset of genes that are informative varies based on the triplet of cell types being considered, establishing lineage relationship between all cell types at once as opposed to three at a time could be challenging. Indeed, our attempt to reconstruct the lineage relationships between all cell types using methods based on PCA failed ( Figure 1D).
We further evaluated the clear minimum pattern in 100 triplets in which there was no clear relation between the cell types ( Figure 1-source data 1). We found that while there were a substantial number of genes exhibiting a clear minimum in one of the three cell types in unrelated triplets, their minima were evenly distributed amongst the cell types. To quantify that the minima were evenly distributed, we counted the fraction of genes f i which reached a clear minimum in cell type i = A, B, C (for A, B, and C unrelated cell types), and for each triplet, we computed the entropy We compared the distribution of the entropy and the number of genes showing a minimum in any triplet for unrelated and related triplets ( Figure 1-figure supplement 3A). The unrelated triplets have higher entropy and typically more genes with a minimum level. This suggested that unrelated triplets show a distinct pattern from the related triplets.

Using patterns to infer lineages
Together, these observations suggested a strategy for inferring the lineage topology between three cell types: each gene showing the clear minimum pattern with a minimum expression level in a particular cell type increases to the probability that this cell type is not the root of the topology ( Figure 2A). We next developed a statistical machinery to systematically detect this pattern in gene expression data and to use the resulting sparse subset of genes to infer lineage relationships between three cell types at a time. We then used the inferred relationships between all sets of three cell types as constraints to determine the full developmental lineage tree. The classifications of genes as transition or marker genes in Figure 1 were based on p<0.005 in a two-sample t-test. To implement such classifications probabilistically without arbitrary cutoffs we developed a statistical framework to infer the lineage relationships between each set of three cell types A, B, and C and find the key sets of transition genes (those genes that show the clear minimum pattern), given gene expression data g A;B;C i n o in those cell types. We determined the probability of any possible topological relationship between the cell types T ¼ A; B; C; f g referring to cell type A, B, or C being the root of the triplet, or which corresponds to the case where the data does not suggest any lineage relationship between the three cell types because either no significant pattern could be detected, or multiple genes exhibiting the minimum pattern suggested conflicting topologies. Rather than an absolute classification of genes as showing a pattern or not, we calculated the  We derived the probability of a given topology T given the expression data as (Materials and methods): is the odds of gene i being a transition gene and thus having a unique minimum. The term p i T is min j g A;B;C i is the probability that the mean i T of the distribution of the expression levels of gene i in the root cell type T is less than the mean in the other two cell types. The odds implicitly contains the only free parameter in our analysis, the prior odds Þ , which defines the number of genes we expect to show the clear minimum pattern a priori, and functions as a sparsity parameter for the inference. Qualitatively, in the above equation, every gene casts a vote Àp i T is min j g A;B;C i against the cell type T in which its mean expression is minimal being the root. Further, this vote is weighted by the odds O i of gene i being a transition gene. Thus, genes with a clearer minimum pattern get larger votes in determining which cell type is not the root. In practice, these quantities are computed numerically (Materials and methods). We note further that if a substantial number of genes cast votes against each of the cell types, then the probability of the null topology increases. We computed the probability of obtaining the null topology among the 150 related triplets and 100 unrelated triplets from our training set. The distribution of the probability of obtaining the null topology was considerably different between the related triplets and the unrelated triplets, with an AUC of 0.96 (Figure 1-figure supplement 3B-C).

Application to hematopoietic gene expression data
We used our statistical framework to recreate the lineage of early hematopoietic differentiation. We considered 11 early hematopoietic progenitors from the ImmGen Consortium microarray data set (Heng et al., 2008) (Figure 2-source data 1). These cell types and their associated relationships were not included in the data set used earlier to study the correlations of the two patterns and lineage topologies. Several features of the early hematopoietic lineage tree are debated against the topology whose root has the minimum mean among the three cell types, and this vote is weighted by the odds that the gene is a transition gene (Equation 1). Two groups of genes, labeled by their names, have high odds of being transition genes and thus cast a strong vote against CMP or MPP being the root.(C) The computed probability of the topology given gene expression data indicates 0.84 probability that ST is the intermediate type. (D) The plot of the probabilities of genes being transition genes for triplet ST/MPP/CMP, given gene expression data and that the topology is MPP -ST -CMP. The names of the 10 genes with the highest probability of being transition genes are shown. Probabilities are calculated assuming There are two classes of transition genes: one for which gene expression in CMP is greater than expression in MPP (regular font), and another for which gene expression in MPP is greater than expression in CMP (bold font).(E) Plot of the replicates of ST, MPP and CMP in the gene-expression space of the two classes of transition genes (with probability > 0.8). Plotted on each axis is the mean normalized log expression level of the transition genes in the class, each class is denoted in curly brackets by the name of the transition gene with the highest probability.(F) Dot plot for triplet MEP/GMP/FrBC of the cell type that is most likely to have the minimum mean expression as a function of the odds O i of that gene being a transition gene.(G) The computed probability of the topology given gene expression data is the null hypothesis (p=0.99). DOI: 10.7554/eLife.20488.008 The following source data and figure supplements are available for figure 2: To illustrate our method, we first described the analysis of the expression data from two such triplets of cell types: CMP/ST/MPP and MEP/GMP/ FrBC ( Figure 2B-G). We then assembled the triplets to form an undirected lineage tree ( Figure 3; Video 1).
Following Equation 1, each gene votes against the topology whose central node has the minimum expression of that gene among the three cell types, and this vote is weighted by the odds that the gene is a transition gene. To illustrate this for the triplet of cell types CMP, ST and MPP, we plotted the topology each gene voted most against, i.e. the topology T for which is the maximum, versus the odds O i of that gene being a transition gene ( Figure 2B). We find two groups of genes that are much more likely to be transition genes than any of the other genes, with values of O i~1 0 2 compared to 10 0 at most for other genes ( Figure  In fact, the intuition in Figure 2B is borne out in the calculation of p T j g CMP;ST;MPP i n o . Using Equation 1 above and assuming a sparsity parameter of 0.05, we calculate that there is an 84% chance that the topology is ST ( Figure 2C; Figure 2-figure supplement 2B). Although gene Irf8 ( Figure 2B, italic font) is strongly downregulated in ST and is expressed at higher levels in CMP and MPP ( Figure 2C), we note that its signal is overwhelmed by the large number of genes downregulated in either CMP or MPP, illustrating the statistical nature of the framework. For each triplet, we evaluated each gene's probability of being a transition or marker gene (Figure 2-source data 3). Figure 2D shows the names and associated probabilities of the 12 genes most likely to be transition genes for the triplet CMP À ST À MPP. The transition genes fall into two groups, corresponding to the two groups in Figure 2B. One group, which includes genes Gata1, Dach1, and Gata2, has higher expression in CMP than in MPP; the other group, which includes Hlf, Tsc22d1, and Hoxa9, has higher expression in MPP. Although the values of the probabilities of the genes being transition genes vary with the value of the sparsity parameter, the relative order of different genes does not change. The genes identified include many genes previously identified as being important for lineage specification (Crispino, 2005;Gazit et al., 2013;Miyawaki et al., 2015). The transition genes we discovered thus not only have gene expression patterns that reflect the lineage decision but also include functionally important genes.
In addition to the transition genes, we identified marker genes (p a i ¼ 1 j g i f g; T ð Þ>0:8) present only in ST (including Mpl and Rai14, consistent with [Solar et al., 1998]) and then symmetrically downregulated in both CMP and MPP (Figure 2-figure supplement 1C). Marker genes for CMP include Srf, Zeb2, Rbpj and Irf8 (consistent with [Goossens et al., 2011;Kurotaki et al., 2013;Ragu et al., 2010;Robert-Moreno et al., 2005;Tamura et al., 2000]); marker genes for MPP include Satb1, consistent with (Satoh et al., 2013). Although these genes were not used to determine the topology, they are good markers for cell types ST, CMP and MPP.
Plotting the cell types using the mean expression levels of the two transition gene class captures the fork in the gene expression space associated with the cell-fate decision ( Figure 2E). In contrast with the PCA analysis of the cell types (

A lineage tree for early hematopoiesis
We next reconstructed the early hematopoietic lineage tree and identified transition genes involved in the different cell state transitions using all the non-null triplet relationships as constraints on the lineage relationships between all cell types. Out of the 165 possible triplets of hematopoietic progenitors, 144 showed one single non-null topology with probability greater than 0.6 over a range of the prior odds from 10 À6 to 10 2 . We next determined an undirected graph that recapitulates all of the individual triplet topologies (note that we are only inferring triplet topologies and are not inferring directionality).  Video 1. Tree-building process for early hematopoietic lineage Animation of the triplet assembly and pruning process for reconstructing the early hematopoietic lineage. For illustrative purposes, triplets (with p>0.6) are successively selected at random (in practice, the assembly and pruning process was performed on all triplets simultaneously; the resulting tree does not depend on the order in which the triplets are selected). The nodes of the current triplet are highlighted in yellow; if a topology is recognized for the triplet, the root is shown in green and the leaves in yellow, and the triplet edges are shown in magenta. If adding the triplet causes another triplet to be pruned, the soon-tobe-pruned (i.e. offending) edge is highlighted in red. The resulting pruned graph is then shown before adding the next triplet. As more triplets are considered, more edges between nodes are added and then pruned, leading to the final tree. DOI: 10.7554/eLife.20488.018 MPP -CMP (Figure 2-figure supplement 1A, left). According to a model suggested by Adolfsson and colleagues (Adolfsson et al., 2005), ST splits into CMP and MPP (Figure 2-figure supplement 1A, right), and the topology should be CMP -ST -MPP. We identify both CMP -ST -MPP and CMP -LT -MPP as the correct topologies, lending support to the Adolfsson model. The Adolfsson and traditional models differ in the topology of 9 triplets. The inferred expected topologies of 8 out of these nine triplets support the Adolfsson model, which led to the identification of the final tree ( Figure 2-figure supplement 2).
The lineage tree that we determined is consistent with the Adolfsson model and contains three additional lineage decisions ( Figure 3A-B). First, CMP gives rise to MEP (megakaryocyte/erythroid progenitor) and GMP (granulocyte/macrophage progenitor). Second, MPP gives rise to MLP (multilineage progenitor), which then splits into the GMP and CLP (common lymphoid progenitor) cell types. In the final lineage decision, CLP gives rise to the ETP (pre-T) and FrA (pre-pro-B) cell types.
For each triplet of cell types along the tree, we identified transition and marker classes of genes. Among the 14 triplets that contained only adjacent cell types, we identified on average 24 marker genes per cell type and 25 transition genes (probability threshold of 0.8, prior odds . Many genes we discovered as belonging with high probability to the transition and marker classes of genes at each lineage decision are known in the literature to be functionally important genes, including classic hematopoietic regulators such as Cebpa, Sfpi1, Gata1, Satb1, Irf8 and Ebf1 (see full tables with references in Figure 3C and Figure 3-source data 1). Additionally, the genes identified include many genes successfully used in hematopoietic reprogramming experiments, including Gata2 and Pbx1 ( Figure 3C). Together these observations suggest that the sparse subspace of transition and marker genes identified by our framework not only allows for accurate reconstruction of the lineage hierarchy but also constitutes a set of candidates for relevant biological functions.
As further validation of the inference method, we compared it to the method proposed by Grun et al. on a single-cell intestinal development data set (Grün et al., 2015(Grün et al., , 2016. We inferred lineages between each cell type based on their cluster identifications, excluding clusters with fewer than 10 cells, and constructed an undirected lineage tree by taking triplets with probability > 0.6 and applying the pruning rule (Figure 3-figure supplement 2A). The only disagreement between the two methods is the progression from crypt base columnar cells (C 2 ) to different populations of Goblet cells (C 4 and C 8 ). Grun et al. hypothesize a C 2 -C 8 -C 4 progression, while we infer the triplet C 8 -C 2 -C 4 with p>0.99, suggesting that the progenitor C 2 gives rise to both differentiated Goblet subpopulations. Both lineage trees are supported by the literature (van der Flier and Clevers, 2009). The high probability transition genes included many factors well known for their roles in tissue homeostasis and development (Figure 3-figure supplement 2B), notably Klf4 (Yu et al., 2012), Atoh1 (VanDussen and Samuelson, 2010), Spdef (Noah et al., 2010), Foxa1/Foxa2 (Ye and Kaestner, 2009), and Tcf3 (Merrill et al., 2001).
Inferred lineage tree for human excitatory neuronal progenitors from in vitro single-cell data over 80 days of differentiation The ease with which single-cell transcriptomic data can be generated (Grün et al., 2015;Jaitin et al., 2014;Macosko et al., 2015;Patel et al., 2014;Paul et al., 2015;Treutlein et al., 2014;Zeisel et al., 2015) presents an opportunity to understand the dynamics of the underlying networks that lead individual cells to their final fate. We studied the differentiation of stem cells both into germ layer progenitors (Jang et al., 2017) and into cortical neurons. To study the latter, we analyzed single-cell gene expression data from 2217 cells from an in vitro differentiation protocol for early human neuronal development (Yao et al., 2017). Briefly, human embryonic stem cells (hESCs) were subjected to a SMAD inhibition-based cortical induction phase, a progenitor expansion phase, and a neural differentiation phase. Single cells were sorted at 12, 26, 54, and 80 days into differentiation, and their gene expression was profiled using the SMART-Seq2 technique (Picelli et al., 2013). In the initial clustering of the single-cell data, dimensionality reduction by PCA (into 15 above-noise components) followed by t-SNE (Van Der Maaten and Hinton, 2008;Satija et al., 2015) showed separation by day and SOX2 expression ( Figure 4A). However, the number of predicted clusters varied depending on the perplexity parameter (Figure 4-figure supplement 1A-B). In addition, no clear lineage or distance relationship among the putative types is immediately apparent from this clustering. Analysis of this data with other recent methods such as Monocle and     Analyzing data from single-cell profiling presents an additional challenge relative to data from population-level profiling because cell types are not previously known and must be inferred from the data. Computationally, it is necessary define a measure of similarity in the gene expression profiles of individual cells so as to be able to cluster them and define cell states. Here again, it is necessary to identify the correct gene subspace to use for clustering. Clustering and determining lineage have typically been performed sequentially and treated as independent problems Trapnell, 2015;Trapnell et al., 2014). However, we found that it is informative to solve both these problems simultaneously.
In our framework, the relevant feature set for clustering is the set of marker and transition transcription factors from the triplets with non-null topologies. But determination of these sets of genes and of the transition topologies depends itself on knowledge of the cluster identities. Following previous work on sparse clustering (Witten and Tibshirani, 2010), we simultaneously determined . . . ; c n f g, sets of transitions T f g, marker genes (a i ¼ 1) and transition genes (b i ¼ 1), given single-cell gene expression data g i f g. Starting from a seed clustering scheme C 0 f g; iterative maximization of the conditional probabilities p T f g; a i f g; b i f gj g i f g; C f g ð Þ and p C f gj g i f g; T f g; a i f g; b i f g ð Þ converges to most likely set C f g; T f g; a i f g; b i f g ð Þ (C) Cell-cell covariance matrix between cells using only the associated high probability marker and transition genes show the final cluster assignments c 0 ; c 2 and c 3 (right) in contrast to using all transcription factors (left). (D) Selected high probability triplets of clusters plotted in the axes defined by two sets of transition gene classes for each triplet. c 1 À c 0 À c 2 (top right, p T ¼ c 0 j g c0;c1;c2 i È É À Á > 0:99), plotted in transition gene class {CEBPG} also including POU3F1, POU3F2, NR2F1, NR2F2, ARX, LIN28A, TOX3, ZBTB20, PROX1 and SOX15, and class {DMRTA1} also including HES1, HES5, FOXG1, PAX6, HMGA2, SOX2, SOX3, SOX9, SOX6, SP8, OTX2, TGIF, ID4, TCF7L2, and TCFL1. c 2 À c 0 À c 3 (top left, p T ¼ c 0 j g c0;c2;c3 i È É À Á ¼ 0:96), plotted in transition gene class {LHX2} also including FEZF2, FOXG1, HMGA1, SP8, OTX1, SOX11, GLI3, SIX3, ETV5, and class {POU3F2} also including GTF2I, HIF1A, ID1, ID3, PROX1, SALL1, SOX21, TCF12, TRPS1, ZHX2.c 1 À c 0 À c 3 (bottom left, p T ¼ c 0 j g c0;c1;c3 i È É À Á >0:99), plotted in transition gene class {FOXO1} also including HMGA2, PAX6, and SOX2, and class {LHX2} also including DMRTA2, HMGA1, ARX, LIN28A, OTX2, LITAF, NANOG, POU3F1, SOX15. c 6 À c 5 À c 7 (bottom right, p T ¼ c 5 j g c5;c6;c7 optimal clusters, lineage topology and sets of transition and marker genes by iteratively selecting transition and marker genes and clustering the data using this set of features ( Figure 4B).
In order to utilize information about developmental distances in clustering, we iteratively maximized a joint distribution, P T; C f g m j g i f g ð Þ, over the developmental tree, T; and the set of clusters of cells, C f g m ¼ c m 1 ; c m 2 ; . . . ; c m K È É ( Figure 4B). Starting with a clustering C f g m ; we first inferred the set of genes a i ¼ 1 f g and b i ¼ 1 f g which were identified in high probability triplets, p Tj g i f g; C f g m ð Þ>0:6, and we then re-clustered in this new subspace to obtain the clusters for the next iteration C f g mþ1 . Initial analysis based on the gap statistic (Tibshirani et al., 2001) suggested that the single-cell gene expression profiles clustered into 20 clusters of cell types. We chose a seed C f g 0 for the iterated clustering procedure by intentionally over-clustering the data into 40 clusters using spectral K-medoids. By being overly discriminative in our initial clustering, we ensured that all genes with differential regulation would be classified as either marker genes or transition genes, and would be preserved in later clustering iterations. We iterated the clustering-inference procedure until the dimension of the re-clustering subspace changed by less than 10% of the total transcription factor space. In this resulting subspace of 469 genes, we finally clustered the cells into the final configura- . . . ; c 7 f g using K-medoids ( Figure 4C, Figure 4-source data 1) where the number of clusters K = 8 was chosen based on the gap statistic (Tibshirani et al., 2001).
In addition to interpreting these individual transition genes defining the major branch splits, we correlated the expression over all predicted transition and marker genes in neuronal clusters to in vivo developmental human data (Miller et al., 2014). The in vivo data comprise a representative range of microarray data sampled from different parts of the developing brain at post-conception week 15, including forebrain proliferative regions, midbrain, and hindbrain. The differences in data acquisition methods (RNAseq vs. microarray, single-cell vs heterogeneous populations) resulted in relatively low correlations overall, but there are clear associations between individual clusters and specific brain regions ( Figure 4E). Specifically, c 1 , maps to the ganglionic eminences, suggesting an interneuron identity, whereas c 5 ; c 6 and c 7 show better mapping to mid-and hindbrain structures, and c 4 appears to be more closely related to neocortex. Overall, this global comparison, combined with the identification of genes with known regional expression, suggests that the inferred clusters from the in vitro data capture the diversity of differentiation into the early stages of the major neuronal lineages ( Figure 4F). These lineage predictions based on our analysis techniques were verified experimentally using viral barcoding in a separate work (Yao et al., 2017).
To estimate the sparsity of the underlying network and to find a minimal subset of genes through which lineage could be inferred, we replicated the analysis while only considering a limited set of genes per triplet. We assembled a collection of 20 triplets with maximal leaf-to-leaf distance of 4 nodes, and non-null inferred topology. For each triplet, we ranked genes based on their odds of being a transition gene, O i ; agnostic of the true topology of the triplet. We then replicated the inference process using only the N genes with the greatest odds ( Figure 4-figure supplement 1E). We found that with as few as 4 genes per triplet, the correct lineage topology could be inferred for all of the triplets. Genes with greatest odds comprise a restricted subset of genes for further experimental investigation. Further, these findings suggest that the dynamics of expression of just four specific genes are sufficient to monitor a particular lineage decision in single cells.

Discussion
Finding an informative subspace of variables for data analysis is a general problem in machine learning, both for regression and clustering (Tibshirani, 1996;Witten and Tibshirani, 2010); the innovation in this paper is to use a statistical pattern learned from known biology to inform this subspace search. The approach we take here is complementary to methods that project expression variability onto coordinates of PCA, ICA or t-SNE maps (Marco et al., 2014;Satija et al., 2015;Trapnell et al., 2014), which are combinations of all variables. Searching for sparse representation of the dynamics has the advantage of providing interpretability and experimental direction (McGibbon and Pande, 2017). Following the dynamics of this small set of high-probability transition genes via fluorescent tagging could allow for the tracking of lineage decisions of individual cells in real time. Further, these genes provide a list of candidates for drivers of fate decisions, and hence a set of experimental hypotheses.
Not all genes give us equal information about the dynamics of differentiation. We discovered that genes showing the clear minimum pattern are most predictive of the sequence of lineage transitions during development. Although our pattern discovery and subsequent lineage reconstruction does not assume any functional role for the clear minimum pattern, we note that this pattern is shown by genes known to be regulators of development during hematopoiesis. The same pattern observed in many differentiating systems (Graf and Enver, 2009;Qi et al., 2013;Thomson et al., 2011;Zhang et al., 1999), and is consistent with mutual inhibition. Mutual inhibition, in turn, is hypothesized to play an important role in maintaining multi-stable systems and in mediating transitions between different stable states of multi-stable systems (Ferrell, 2012).
Discovery of sparse representations of the cell states and variability between them demonstrates the efficacy of low dimensional descriptions of the system. Understanding the dynamics and transitions of complex physical systems composed of a large number of variables has been driven by the discovery of low dimensional order parameters (Anderson, 1978;Landau and Lifshitz, 1951). As opposed to measuring and modeling states as high dimensional objects in their native representation, order parameters provide low dimensional descriptions of the states and dynamics, which has proven crucial in developing both qualitative and quantitative models. Finding small subsets of genes which captures the lineage transitions in cells analogously provides a low dimensional subspace that captures the dynamics in genetic networks and can be useful for modeling (Jang et al., 2017). An accompanying paper allows us to exploit this idea to extract mathematical models for the underlying molecular circuits from single cell gene expression data obtained during germ layer differentiation (Jang et al., 2017).

In vitro neuronal differentiation
Single-cell transcriptomic data from the in vitro neural differentiation procedure was obtained as described in Yao et al. (2017) hESCs were dissociated with Accutase and plated on Matrigel-coated 24-well plates at 2.5 Â 105 cells/cm2 in DMEM/F12 (#11330-032), 1 Â N2, 1 Â B27 without vitamin A, 2 mM Glutamax, 100 mM non-essential amino acids, 0.5 mg/mL BSA, 1X Pen-Strep, and 100 mM 2-mercaptoethanol (referred to as basal medium; all from Thermo Fisher, Waltham, MA) with 20 ng/mL FGF2 (Thermo Fisher) and 2 mM thiazovivin. Cortical induction was initiated by changing to the basal medium with 5 mM SB431542 (StemRD, Burlingame, CA), 50 nM LDN193189 (Reagents Direct, Encinitas, CA) and 1 mM cyclopamine (Stemgent, Lexington, MA) (referred to as NIM). NIM was changed daily for 11 days. On day 12, cells were dissociated and seeded on Matrigelcoated 24-well plates at 5 Â 105/cm2 in basal medium with 20 ng/mL FGF2 and 2 mM thiazovivin. Progenitor expansion was initiated on D13 by changing to serum-free human neural stem cell culture medium (NSCM, #A10509-01 from Thermo Fisher) containing 20 ng/mL FGF2 and 20 ng/mL EGF. NSCM was changed daily for 6 days. Cultures were passaged once more on D19 with Accutase and replated at 5 Â 105 cells/cm2. On D26, cells were dissociated with Accutase and seeded on 24-well plates sequentially coated with poly-D-lysine (Millipore, Billerica, MA) and laminin (Thermo Fisher) at 1 Â 105 cells/cm2 in basal medium supplemented with 20 ng/mL FGF2 and 2 mM thiazovivin. On D27, medium was changed to a 1:1 mixture of DMEM/F12 and Neurobasal medium (#21103-049) supplemented with 100 mM cAMP (Sigma-Aldrich, St Louis, MO), 10 ng/mL BDNF (R and D Systems, Minneapolis, MN), 10 ng/ mL GDNF (R and D Systems) and 10 ng/mL NT-3 (R and D Systems) (referred to as ND). Cells were maintained in ND medium for four weeks until day 54 with half medium change every other day. Quality of differentiations was routinely assessed by immunostaining at D12 (PAX6 and DCX), at D26 (LHX2, SOX2, EOMES, POU3F2, and TBR1), and at D54 (MAP2 costained with TBR1, CTIP2, SATB2). In addition flow cytometry at D26 (EOMES, SOX2 and PAX6) was performed. Typically, EOMES at day 26 proved the most valuable quality control metric (~10% of cells by both flow cytometry and immunostaining) and predicted failure at D54. Specifically, when EOMES was low (<1% of cells) differentiations failed and were typically dominated by POU3F2+ cell types and/or non-neural 'other' cell types. These failed differentiations were eliminated from further analysis, typically~20% of experiments (5 of 19 experiments in 2016). 50 bp paired-end Smart-Seq2 libraries were prepared from these cells as previously described (Picelli et al., 2013) and mapped as described in Thomsen et al., 2016 andYao et al. (2017). As each cell was profiled independently (without pooling before amplification, as in methods such as Cel-Seq, STRT, or Drop-Seq), we did not observe the batch effects present in pooling-based methods. Although cells were profiled in plates (batches), there was no significant plate-related variation -this Fixed single-cell transcriptomic characteriza is most likely due to amplification being carried out separately in each well of the plate, thereby precluding cross-talk among barcodes, or plate-related differences in amplification. This is clear from the mixing of cells from different plates in the clustering. A complete census is provided (Figure 4source data 4).

Gene expression data
Hematopoietic gene expression data were downloaded from the Immunological Genome Project (Heng et al., 2008); GEO GSE15907) and log-2 transformed. We restricted the genes considered to 1459 transcription factors. Brain development expression data were log-2 transformed, and 1460 transcription factors were individually normalized by dividing by the 90th-percentile expression value. Lists of mouse and human transcription factors are provided (Figure 1-source

Description of algorithm
The algorithm proceeds according to the following steps: 1. Find initial seed clustering configuration C f g 0 using K-medoids where K is chosen to be larger than the cluster number inferred from the Gap statistic (Tibshirani et al., 2001) 2. For all triplets of clusters, find most likely T and a i f g and b i f g given C f g m : 4. Repeat steps 1 and 2 until convergence of C f g: 5. Determine the most likely tree connecting cell clusters, recapitulating high-probability triplet topologies.
Each of these steps is described in the following section.
Bayesian framework for inferring cluster identities, state transitions, and marker and transition genes simultaneously Notation; Bayes' rule Given gene expression data from single cells g i f g, we built a probabilistic framework to simultaneously infer cell cluster identities, C f g c A ; c B ; . . . f g, the sequence of transitions T between these clusters, the key sets of marker genes a i f g that define each cell cluster, and genes b i f g that determine the sequence of transitions between clusters. We maximized the joint probability distribution of these variables given the gene expression data, p T; c A ; c B ; . . . f g; a i f g; b i f g j g i f g ð Þ to determine the maximum likelihood estimates of these parameters.
We first consider how to solve this problem in the case in which there are three cell clusters, and we will later build a tree using all possible combinations of three cell clusters. Let the set of three cell clusters be c A , c B and c C with gene expression data g A;B;C i n o for all genes i ¼ 1 to N ð Þ and all cells. The term g A;B;C i denotes the expression data for just gene i in cells in clusters c A , c B , and c C . The topology T of the relationships between cell clusters c A , c B and c C can take on four possible values: T ¼ A: cell cluster c A is in the middle (either c A is the progenitor of c B and c C , or c A is an intermediate cell type between c B and c C ); T ¼ B: cell cluster c B is in the middle; T ¼ C: cell cluster c C is in the middle; or T ¼ : we cannot determine the topology. Complementarily, for each gene i we define variables a i and b i , where a i ¼ 1 and b i ¼ 0 if the gene is a marker gene, a i ¼ 0 and b i ¼ 1 if the gene is a transition gene, and a i ¼ b i ¼ 0 otherwise. Our task is to determine the probability p T; c A ; c B ; . . . f g; a i f g; b i f g j g A;B;C i n o given gene expression data for all genes g A;B;C i n o .
According to Bayes' rule, p T; c A ; c B ; . . . f g; a i f g; b i f g j g A;B;C i n o is proportional to the probability of the gene expression data given T, C f g; a i f g and b i f g: The denominator of the right hand side of Equation 2 is a normalization constant.
and p C f g ð Þ are respectively the prior probabilities of a i f g and b i f g given T and C f g; the prior probability of T given C f g, and the prior probability of C f g. We assume that in the absence of any expression data, the probability that a gene is a transition or marker gene is independent of the topology and clustering configuration:

Conditional independence
In our model, we assume that knowing the clustering configuration C f g, the topology T and whether or not a gene is a marker or transition gene is sufficient to determine the probability distribution for its expression levels in each of the cell clusters. Therefore, the gene expression patterns of different genes are conditionally independent given the topology, clustering and gene type: Thus, Equation 2 becomes We maximize the evaluated p T; c A ; c B ; . . . f g; a i f g; b i f g j g A;B;C i n o with respect to T, C f g and each of the a i and b i to obtain the most likely relationships between cell types c A , c B and c C , as well as the genes most likely to be marker and transition genes.  in the three cell types A, B and C has the smallest mean value in either B or C but not in A (Figure 2A). If the distribution of the expression of gene i in cell type A is D A g A i j i A ; s i A ; C f g À Á (we assume a log-normal distribution) with a mean i A and standard deviation s i A , with analogous expressions for cell types B and C, then our model defining p g A;B;C

Expression for p g
where i A , i B and i C are the mean values of the expression levels of g i in cell types A, B and C. Thus, The terms in Equation 5 can be calculated by integrating over the prior probability distribution of the means i A , i B and i C and standard deviations s i A , s i B and s i C , with the conditions on the means constraining the domains of integration: In addition to topologies A, B and C, we consider a null hypothesis in which transition genes have differential expression levels between states, but these levels are not correlated with any particular topology of states. This corresponds to having gene expression levels from cell-types A, B and C coming from three distributions with no restrictions on the relative order of the three means: Note that the probability of the data given the null hypothesis is the average of the probabilities of the data given the non-null hypotheses: Note that pðg A;B;C i jT; C f g; a i ¼ 0; b i ¼ 1Þ depends on both T and C f g.
Expression for p g A;B;C i j T; C f g; a i ¼ 1; b i ¼ 0 (marker genes) Our model for marker genes assumes that the probability distribution for the expression level of such genes, p g A;B;C i j T; C f g; a i ¼ 1; b i ¼ 0 to be independent of T and to be generated from distributions with two cell-types having a low value and the third a high value (for example, for cell-types A and B and D C g C i j i C ; s i c À Á for cell-type C, with the constraint i AB < i C Þ: Note that p g A;B;C i j T; C f g; a i ¼ 1; Our model for genes that are neither marker nor transition genes is that the expression levels of such genes, p g A;B;C

Numerical integration
Each of the probabilities on the right hand side of Equation 4 is evaluated numerically as above. We assume the distribution of the expression of gene i in cluster Given m log2-transformed replicate measurements g A i of gene expression of gene i in cells belonging to cluster c A , the probability of the data assuming mean i A and standard deviation s i A is: Distributions D B ; D C ; D AB ; D AC ; D BC and D ABC are defined analogously. We take the a priori probability distribution of i and s i , p i ; s i ð Þ as uniform over a certain range of means and standard deviations estimated from the data. For the log2-transformed hematopoietic gene expression data, we take 2< i <14 and 0<s i <0:75.
We take the prior probabilities for the distributions in different cell types to be independent: The constraints on the order of the means are enforced by the domain of integration, and the prior must be properly normalized over this domain. For example, in Equation 6, Integrals were evaluated numerically in MATLAB using trapezoidal integration with step-sizes d = 0.05 and ds = 0.01.
The choice of a log-normal distribution could potentially confound the results, particularly for single-cell RNA-Seq data with significant zero-inflation. This is a potential area for improvement in our algorithm, but the method can easily be adapted to model different distributions of RNA expression, such as gamma distributions (Shahrezaei and Swain, 2008;Wills et al., 2013) or beta-Poisson distributions (Delmans and Hemberg, 2016;Vu et al., 2016). In either case, the probability of the data given different topologies would be computed by numeric integration over the parameters of the distribution, for example, a and b for the Gamma distribution, by replacing the log-normal distributions in Equations 7, 9 and 10 with ones from the appropriate model.
The choice of the right parametric form for single-cell RNA expression is still an area of active research. Our choice of log-normal distributions assumes that higher order moments of the distributions (beyond standard deviation) ought to have a minimal contribution to the predictions, but we have not tested this extensively.
Although the default prior was the uniform prior, we also implemented an empirical prior p ; s ð Þ by estimating the empirical distribution over all the Immgen cell types, using the kernel density estimation code provided in (Botev et al., 2010). The resulting hematopoietic lineage tree was identical. Using the kernel-density-estimated empirical prior may provide more stability in future analyses.
Probability of topology given gene expression and cluster identities We can derive the probability of the topology given the gene expression data and cluster identities p T j g A;B;C i n o ; C f g by summing over all the a i f g and b i f g to find the probability of the data given topology p g A;B;C i n o j T; C f g : where the probability of gene expression data for gene i given topology p g A;B;C i j T; C f g is obtained by summing p g A;B;C i j T; C f g; a i ; b i over a i and b i : We consider the odds of a gene being a transition gene Þ as a sparsity parameter, which we vary. We take p a i ¼ 0 The probability of topology T given data is proportional to the probability of the data given topology T (using Bayes' rule): where Therefore, using Equation 14, we obtain the following expression for p T j g A;B;C i n o ; C f g :

Rewriting Equation 21 in terms of negative votes
Let us denote the probability of gene expression data for gene i given that cell cluster has the distribution with minimum mean expression as p g A;B;C i j i is min; C f g : For example, where p a i ; b i j T; C f g; g A;B;C i is the probability that gene i is a marker or transition gene given its gene expression, the clustering, and that the topology is T.

Choice of prior odds does not affect the most likely topology
The only free parameter in our calculation above is the prior odds of gene i being a transition gene, then the null hypothesis dominates: if all genes are transition genes, then there will be negative votes against all topologies. We computed the behavior of p T j g A;B;C i n o between these two limits to determine the sensitivity of our answer to Þ between 10 À2 and 10, whereas for triplet MEP/GMP/FrBC there is no value of the prior odds that strongly favors a non-null topology. For most triplets, the most likely topology does not depend on the choice of prior odds; when building lineage trees, we ignore those triplets where different choices of prior odds lead to different most-likely topologies, i.e. there is more than one non-null topology that reaches probability 0.6 over the range of prior odds.

Determination of lineage tree from triplet topologies Selection of triplets
In order to build lineage trees from the topologies we determine for each cell type, we select the triplets for which our determination of the topology is most robust. There is one free parameter in our model: the prior odds for a gene to be a transition gene in the absence of gene expression data, For each triplet, we vary this parameter between 10 À6 and 10 2 and calculate the probability of the topology given gene expression data p T j g A;B;C i n o as a function of the prior odds. We want to consider only triplets which showed a single dominant topology. We exclude triplets which show a weak probability for a particular topology or ones which depend on a particular choice of prior odds. We also do not consider triplets which show a strong probability for two different topologies, depending on the choice of prior odds.
There were 14 such triplets among the 165 hematopoietic triplets. Mathematically, these cases come up when the genes that are most likely to show the clear minimum pattern (furthest on the right in a 'dot plot') suggest one topology, but if one used a more permissive value of the sparsity parameter, a different topology wins out. One of the cell types might have a small number of genes with very high odds, but then fewer genes with moderately high odds compared to the other cell types. We did not notice a clear pattern in the identity of the triplets exhibiting this behavior, but 9 of the 14 were triplets of length five or greater in the Adolfsson model. One of the triplets with this behavior was the MLP/CMP/GMP triplet, and the dominant topology was either MLP (at low prior odds) or CMP (at higher prior odds). Interestingly, both cell types are progenitors to GMP in the Adolfsson model. We consider triplets for which only one non-null topology has probability p T j g A;B;C i n o greater than 0.6. The probabilities of the different topologies for each triplet in the hematopoietic tree are shown in Figure 2-source data 1 and for the triplets in the cortical development tree in Figure 4source data 2.

Pruning rule
We assemble the triplets with known topology into an undirected graph. Since we determined topologies by considering cell types three at a time, we obtain topological relationships involving both cell types that are nearest neighbors and cell types that are more distantly related. In order to reconstruct the tree, we must determine which cell types are nearest neighbors and which ones are separated by one or more intermediate cell types. The set of inferred topologies allows us to determine which cell types are separated by intermediates. For every pair of cell types, we ask whether any of the inferred topologies features an intermediate between the two cell types. If such a topology has been inferred, we consider that the two cell types are not nearest neighbors, and that at least one other cell type is an intermediate. For example, we can ignore triplet CMP -LT -MPP because triplet LT -ST -CMP testifies that there exists an intermediate between LT and CMP, and triplet LT -ST -MPP testifies that there exists an intermediate between LT and MPP (Figure 3-figure supplement 1).
Note that this pruning rule does not assume the absence of loops. The lineage tree we infer for the hematopoietic progenitors contains a loop that includes ST to CMP to GMP on one side and ST to MPP to MLP to GMP on the other side. The loop involves triplets CMP -ST -MPP, ST -CMP -GMP, ST -MPP -MLP and MPP -MLP -GMP (we cannot determine the topology of triplet CMP/ MLP/GMP). None of the triplets shows a topology that would allow us to break up the loop.

Distinguishing between two models of hematopoiesis
The topologies we infer support the model from Adolfsson et al. (2005), in which CMP splits from ST-HSC. In particular, there are several triplets that can distinguish the Adolfsson model from the traditional picture, and they support the Adolfsson model.

Stability analysis
The inference algorithm depends on several parameters and priors. We performed a stability analysis for both the microarray hematopoietic data and the human brain single cell data to determine the parameter ranges for which the inferred lineage tree was unchanged.

Hematopoiesis
. For the prior probability, given that a gene is not a transition gene, that it is an irrelevant gene, our default value was p a i ¼ 0 j b i ¼ 0 ð Þ =0.5, but the tree was unchanged for values of p a i ¼ 0 j b i ¼ 0 ð Þ between 0.25 and 0.65.
. For the prior probabilities of different topologies, our default value was p ð Þ ¼ p A ð Þ ¼ p B ð Þ ¼ p C ð Þ ¼ 0:25. We varied p ð Þ while keeping the prior probabilities of the non-null topologies equal: p T 6 ¼ ð Þ ¼ 1 3 1 À p ð Þ ð Þ. The tree was unchanged for values of p ð Þ between 0.1 and 0.35.
. We used a threshold of 0.6 to consider a triplet significant for the tree-building step. The tree was unchanged for thresholds between 0.5 and 0.65.
. A key input parameter into the algorithm is the expected prior distribution of means and standard deviations p ; s ð Þ used in the numeric integration. Our default prior was a uniform prior for both means and standard deviations over reasonable ranges of these parameters. We also implemented an empirical prior p ; s ð Þ by estimating the empirical distribution over all the Immgen cell types, using kernel density estimation (Botev et al., 2010). The resulting hematopoietic lineage tree was identical. Using the kernel-densityestimated empirical prior may provide more stability in future analyses. 2. Brain Development . For the prior probability, given that a gene is not a transition gene, that it is an irrelevant gene, our default value was p a i ¼ 0 j b i ¼ 0 ð Þ = 0.5, but the tree was unchanged for values of p a i ¼ 0 j b i ¼ 0 ð Þ between 0.02 and 0.97.
. For the prior probabilities of different topologies, our default value was p ð Þ ¼ p A ð Þ ¼ p B ð Þ ¼ p C ð Þ ¼ 0:25. We varied p ð Þ while keeping the prior probabilities of the non-null topologies equal: p T 6 ¼ ð Þ ¼ 1 3 1 À p ð Þ ð Þ. The tree was unchanged for values of p ð Þ between 0.06 and 0.96.
. We used a threshold of 0.6 to consider a triplet significant for the tree-building step. The tree was unchanged for thresholds between 0.5 and 0.8.
. For the single cell data, we also changed the number of initial seed clusters for the iterative algorithm within a range from 30 to 45. In each case, the tree remained the same, and never more than 23% of the cells were clustered into a different cell type.