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.

Figure 1 with 3 supplements see all

Download asset Open asset

Figure 1—source data 1

Triplets used for pattern detection.

Sheet 1: The 150 related triplets used for pattern detection are listed, with the root cell type of the triplet in the first column and the two leaves in second and third columns. Cell type names follow (Heng et al., 2008). For each triplet are listed the number of genes with a clear minimum in each cell type and a clear maximum in each cell type (t-test p-value < 0.005). Also listed for each triplet are the length of the triplet, whether the triplet is a lineage progression or a cell fate decision, whether the triplet contains a terminal node, and the entropy in the clear minimum pattern and clear maximum pattern. Sheet 2: The 100 unrelated triplets are listed; cell type names follow (Heng et al., 2008). For each triplet are listed the number of genes with a clear minimum in each cell type and a clear maximum in each cell type (t-test p-value < 0.005) and the entropy in the clear minimum pattern and clear maximum pattern.

https://doi.org/10.7554/eLife.20488.003

Download elife-20488-fig1-data1-v1.xlsx

Figure 1—source data 2

List of mouse transcription factors.

The list of human transcription factors used to cluster and infer lineages for the hematopoietic and intestinal development trees. The list was adapted from (Thomson et al., 2011).

https://doi.org/10.7554/eLife.20488.004

Download elife-20488-fig1-data2-v1.xls

We identified 150 triplets of cell types with experimentally verified lineage relationships from B- and 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 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; Figure 1—figure supplement 1A); the clear maximum pattern: the gene has a clear maximum level in one of the cell types (p<0.005 in both sample t-tests; Figure 1—figure supplement 1B). Note that both patterns can be satisfied simultaneously if the distribution of expression levels in the three cell types are all well-separated with a clear maximum and minimum. We tested if either of the two patterns correlated with the lineage topologies between the three cell types with known lineage relationships (Figure 1C, Figure 1—figure supplement 1C).

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 (Figure 1—figure supplement 1A; Figure 1—figure supplement 2C). Thus the expression pattern of the total set of clear minimum genes correlated with the topology. Since genes showing the clear minimum pattern correlated with lineage relationships (Figure 1C) between cell states both in the case of branches and linear sequences of cell state transitions, we refer to them as transition genes.

We further verified that the clear minimum pattern could be observed (a) in the set of all 25,194 genes (Figure 1—figure supplement 2A), (b) using FDR-adjusted p-values (Benjamini and Hochberg, 1995) (Figure 1—figure supplement 2B), (c) in triplets of different lengths along the lineage tree (Figure 1—figure supplement 2D), (d) in triplets both containing only internal nodes and including terminal nodes (Figure 1—figure supplement 2E), (e) and in both lineage progression and cell fate decision triplets (Figure 1—figure supplement 2F).

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 ${\displaystyle S=-{\sum }_{i=A,B,C}{f}_i\log \left(f_i\right)}$. 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.

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.

Figure 2 with 2 supplements see all

Download asset Open asset

Figure 2—source data 1

Early hematopoietic cell types considered.

Listed for each early hematopoietic progenitor whose name is abbreviated in this paper are the Immunological Genome Project descriptor for the cell type, its common name and phenotype, its age and location, and the number of replicates in the data set.

https://doi.org/10.7554/eLife.20488.009

Download elife-20488-fig2-data1-v1.docx

Figure 2—source data 2

Probabilities of topologies for triplets of hematopoietic cell types.

Listed are, for each triplet of cell types, the probabilities of the four topologies for prior odds $\frac{p\left({\beta }_i=1\right)}{p\left({\beta }_i=0\right)}=0.05$; the number of topologies that reach probability ${\displaystyle p\left(T|\left\{g_i^{A,B,C}\right\}\right)>0.6}$ for some value of $\frac{p\left({\beta }_i=1\right)}{p\left({\beta }_i=0\right)}$ between 10⁻⁶ and 10²; the non-null topology that has the highest probability ${\displaystyle p\left(T|\left\{g_i^{A,B,C}\right\}\right)}$ over the range of prior odds (if the null topology is the most likely topology for the entire range of prior odds, the topology is marked ‘null’); and the value of highest probability ${\displaystyle p\left(T|\left\{g_i^{A,B,C}\right\}\right)}$ over the range of prior odds; the correct topology and triplet length in the traditional model; and the correct topology and triplet length in the Adolfsson model.

https://doi.org/10.7554/eLife.20488.010

Download elife-20488-fig2-data2-v1.xlsx

Figure 2—source data 3

Probabilities of transition and marker genes for the hematopoietic lineage tree.

Listed are, for crucial triplets along the lineage tree, the 200 genes with the highest probabilities of belonging to the two transition gene classes and the root marker class, and their associated probabilities.

https://doi.org/10.7554/eLife.20488.011

Download elife-20488-fig2-data3-v1.xlsx

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 $\left\{g_i^{A,B,C}\right\}$ in those cell types. We determined the probability of any possible topological relationship between the cell types $T= \{A, B, C, \varnothing \}$ referring to cell type A, B, or C being the root of the triplet, or $\varnothing$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 probability ${\displaystyle p\left({\alpha }_i=1|g_i^{A,B,C}\right)}$ of each gene $i$ being a marker gene (denoted by ${\alpha }_i=1$), i.e. gene $i$ showing the clear maximum pattern, and the probability ${\displaystyle p\left({\beta }_i=1|g_i^{A,B,C}\right)}$ of it being a transition gene denoted by ${\beta }_i=1$, i.e. gene $i$ showing the clear minimum pattern (Materials and methods).

We derived the probability of a given topology ${\displaystyle T}$ given the expression data as (Materials and methods):

where, ${\displaystyle {\mathcal{𝒪}}_i=p\left({\beta }_i=1|g_i^{A,B,C}\right)/p\left({\beta }_i=0|g_i^{A,B,C}\right)}$ is the odds of gene $i$ being a transition gene and thus having a unique minimum. The term ${\displaystyle p\left({\mu }_T^i\mathrm{i}\mathrm{s}\mathrm{m}\mathrm{i}\mathrm{n}|g_i^{A,B,C}\right)}$ is the probability that the mean ${\mu }_T^i$ 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 $\frac{p\left({\beta }_i=1\right)}{p\left({\beta }_i=0\right)}$, 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 ${\displaystyle -p\left({\mu }_T^i\mathrm{i}\mathrm{s}\mathrm{m}\mathrm{i}\mathrm{n}|g_i^{A,B,C}\right)}$ against the cell type $T$ in which its mean expression is minimal being the root. Further, this vote is weighted by the odds ${\displaystyle {\mathcal{𝒪}}_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 $\varnothing$ 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).

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 (Adolfsson et al., 2005; Iwasaki and Akashi, 2007) (Figure 2—figure supplement 2A). Given only the gene expression data for these different subpopulations of cells, we determined the lineage relationships and the key factors associated with each lineage decision. We calculated the probabilities of topology and marker and transition genes for the $\left(\begin{array}{c}11\\ 3\end{array}\right)=165$ possible triplets of cell types using our statistical framework (Figure 2—source data 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).

Figure 3 with 2 supplements see all

Download asset Open asset

Figure 3—source data 1

Marker genes for early hematopoiesis.

Table with selected marker genes for early hematopoietic cell types, along with references to published validations of their functional role. Genes known to be effective for reprogramming are shown in bold.

https://doi.org/10.7554/eLife.20488.015

Download elife-20488-fig3-data1-v1.docx

Video 1

Download asset

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 ${\displaystyle p\left({\mu }_T^i\mathrm{i}\mathrm{s}\mathrm{m}\mathrm{i}\mathrm{n}|g_i^{CMP,ST,MPP}\right)}$ is the maximum, versus the odds ${\displaystyle {\mathcal{𝒪}}_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\sim {10}^2$ compared to ${10}^0$ at most for other genes (Figure 2B, regular and bold fonts). These two groups of genes have a large value for either ${\displaystyle p\left({\mu }_{CMP}^i\mathrm{i}\mathrm{s}\mathrm{m}\mathrm{i}\mathrm{n}|g_i^{CMP,ST,MPP}\right)}$ or ${\displaystyle p\left({\mu }_{MPP}^i\mathrm{i}\mathrm{s}\mathrm{m}\mathrm{i}\mathrm{n}|g_i^{CMP,ST,MPP}\right)}$ and thus vote against ${\displaystyle T=\mathit{C}\mathit{M}\mathit{P}}$ (cell type CMP is the intermediate) or against ${\displaystyle T=\mathit{M}\mathit{P}\mathit{P}}$ (cell type MPP is the intermediate). Together these genes that have a high odds of being transition genes appear to most support topology ${\displaystyle T=\mathit{S}\mathit{T}\equiv CMP-\mathit{S}\mathit{T}-MPP}$.

In fact, the intuition in Figure 2B is borne out in the calculation of ${\displaystyle p\left(T|\left\{g_i^{CMP,ST,MPP}\right\}\right)}$. Using Equation 1 above and assuming a sparsity parameter of 0.05, we calculate that there is an 84% chance that the topology is ${\displaystyle \mathit{S}\mathit{T}}$ (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 ${\displaystyle \mathrm{C}\mathrm{M}\mathrm{P}-\mathbf{S}\mathbf{T}-\mathrm{M}\mathrm{P}\mathrm{P}}$. 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 (${\displaystyle p\left({\alpha }_i=1|\left\{g_i\right\},T\right)>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 (Figure 2—figure supplement 1E), in which MPP appears to be an intermediate between the hematopoietic stem cell types (LT and ST) and CMP, the projection of the cell types onto the transition-gene subspace shows that ST splits into CMP and MPP.

In contrast to the case of the triplet of cell types CMP/ST/MPP, for triplet MEP/GMP/FrBC, the distributions of genes supporting different topologies are similar (Figure 2F). Thus the most likely topology calculated using Equation 1 is the null hypothesis (99%), which is that transition genes, if they exist, do not have patterns that depend on the cellular topology (Figure 2G, Figure 2—figure supplement 1D), in which case there is insufficient evidence to classify the triplet according to a particular non-null topology. The distribution of the maximal probability of non-null topologies in the different triplets is heavily concentrated near 1, allowing for clear separation between null and non-null triplets (Figure 2—figure supplement 1E). Null topologies were identified by the algorithm for triplets with cells that are from three terminal nodes (for example, the triplet MEP/GMP/FrBC) or from triplets that contain very distantly related triplets (for example, LT/CLP/ETP) (Figure 2—source data 2).

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⁻⁶ to 10². 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). For example, although triplet CMP/LT/MPP has topology CMP – LT – MPP (Figure 3—figure supplement 1A), we could determine that LT cannot be the direct progenitor of CMP or MPP, because ST is an intermediate between LT and both cell types (Figure 3—figure supplement 1B–C). We could thus 'prune' this triplet when inferring the full graph (Figure 3—figure supplement 1D–E). A visualization of the pruning process is shown in Video 1, where successive triplets are added to the graph, creating new edges and pruning others, leading to the final undirected tree. In practice, though, the pruning process was performed on all triplets simultaneously, not in succession.

The ST/CMP/MPP triplet (Figure 2B–E) immediately distinguishes between two competing models regarding the hierarchy of early hematopoietic progenitors. According to the traditional picture (Iwasaki and Akashi, 2007), MPP is the progenitor of CMP, and ST is the progenitor of MPP – therefore MPP should be an intermediate between ST and CMP and the topology of triplet ST/MPP/CMP should be ST – 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 ${\displaystyle p\left({\beta }_i=1\right)/p\left({\beta }_i=0\right)=0.05}$). 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, 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₂) to different populations of Goblet cells (C₄ and C₈). Grun et al. hypothesize a C₂ – C₈ – C₄ progression, while we infer the triplet C₈ – C₂ – C₄ with p>0.99, suggesting that the progenitor C₂ 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).

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 Monocle2 (Trapnell et al., 2014), TSCAN (Ji and Ji, 2016), and StemID (Grün et al., 2016) did not clearly reconstruct lineage or infer key genes regulating transitions (Figure 4A – Bottom, Figure 4—figure supplement 2). Monocle2 (Figure 4—figure supplement 2A) produces a tree with complex branching, but SOX2+ progenitors and DCX+ differentiated neurons do not clearly separate.

Figure 4 with 2 supplements see all

Download asset Open asset

Figure 4—source data 1

Final cluster identities of single cells from in vitro cortical differentiation.

Cells were assigned to a cluster based on an iterative clustering procedure. Their final cluster assignments are shown here. Cell types 8–11 contained bulk brain control cells, and non-neuronal cells which were excluded from the analysis.

https://doi.org/10.7554/eLife.20488.020

Download elife-20488-fig4-data1-v1.xlsx

Figure 4—source data 2

Probabilities of topologies for triplets of single-cell clusters.

Listed are, for each triplet of clusters, the probability of a given topology and the root of the inferred topology. ‘0’ refers to the null topology. Triplets with p>0.6 were used to construct the lineage tree.

https://doi.org/10.7554/eLife.20488.021

Download elife-20488-fig4-data2-v1.xls

Figure 4—source data 3

Probabilities of Transition and Marker Genes for the Human Brain Developmental Lineage Tree.

Listed are, for crucial triplets along the lineage tree, the genes with the p>0.8 of belonging to the two transition gene classes and the root marker class, and their associated probabilities.

https://doi.org/10.7554/eLife.20488.022

Download elife-20488-fig4-data3-v1.xlsx

Figure 4—source data 4

Human Brain Development SmartSeq2 Census.

Further detail can be found Materials and methods and in Yao et al. (2017) (Supplemental Information).

https://doi.org/10.7554/eLife.20488.023

Download elife-20488-fig4-data4-v1.xls

Figure 4—source data 5

List of Human Transcription Factors.

List of human transcription factors used to cluster and infer lineages for the cortical differentiation tree. List was adapted from (Ben-Porath et al., 2008).

https://doi.org/10.7554/eLife.20488.024

Download elife-20488-fig4-data5-v1.xls

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 (Satija et al., 2015; 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 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, ${\displaystyle P\left(T,{\left\{C\right\}}^m|\left\{g_i\right\}\right)}$, over the developmental tree, $T,$ and the set of clusters of cells, ${\displaystyle {\left\{C\right\}}^m=\left\{c_1^m,c_2^m,\dots ,c_K^m\right\}}$ (Figure 4B). Starting with a clustering ${\left\{C\right\}}^m,$ we first inferred the set of genes $\left\{{\alpha }_i=1\right\}$ and ${\displaystyle \left\{{\beta }_i=1\right\}}$ which were identified in high probability triplets, ${\displaystyle p\left(T|\left\{g_i\right\},{\left\{C\right\}}^m\right)>0.6}$, and we then re-clustered in this new subspace to obtain the clusters for the next iteration ${\left\{C\right\}}^{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 ${\left\{C\right\}}^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 configuration ${\displaystyle {\left\{C\right\}}^f=\left\{c_0,c_1,\dots ,c_7\right\}}$ 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).

We inferred 45 high probability triplets (${\displaystyle p\left(T|\left\{g_i\right\}\right)>0.6}$) between the final cell clusters (Figure 4—source data 2). Four such triplets are shown in Figure 4D, plotted in axes defined by transition genes for each triplet (Figure 4—source data 3). Starting with progenitor cell states (SOX2+, Figure 4—figure supplement 1C), we manually appended cell clusters to the tree according to their time information and in agreement with inferred topological restrictions. The first transition involves the production of day 12 neuronal cell type $c_1$ from day 12 progenitor $c_0$, which is observed in the triplets $c_1-c_0-c_2$ and $c_1-c_0-c_3$. The ${\displaystyle c_1-c_0-c_2}$ triplet (${\displaystyle p\left(T=c_0|\left\{g_i^{c_0,c_1,c_2,}\right\}\right)=0.99}$, Figure 4D – top left) is mediated by 47 transition genes between $c_0$ and $c_1$ and 87 between $c_0$ and $c_2$ (${\displaystyle p\left({\beta }_i=1|\left\{g_i\right\},T\right)>0.8}$). Transition genes expressed highly in the $c_1-c_0$ branch include CEBPG, POU3F1, POU3F2, NR2F1, NR2F2, ARX, LIN28A, TOX3, ZBTB20, PROX1, and SOX15 which have been previously implicated in proliferation of forebrain progenitors (Au et al., 2013; Borello et al., 2014; Cimadamore et al., 2013; Dominguez et al., 2013; Yang et al., 2015) and the migratory behaviors of ganglionic eminences (Kanatani et al., 2008; Kessaris et al., 2014; Lodato et al., 2014; Olivetti and Noebels, 2012; Reinchisi et al., 2012). The $c_0-c_2$ transition genes include DMRTA1, HES1, HES5, FOXG1, PAX6, HMGA2, SOX2, SOX3, SOX9, SOX6, SP8, OTX2, TGIF, ID4, SOX3, TCF7L2, and TCFL1 which are known to be expressed in forebrain progenitors of the developing neocortex, thalamus, and hypothalamus (Abraham et al., 2013; Pozniak et al., 2010; Shimojo et al., 2011; Tzeng and de Vellis, 1998; Wang et al., 2006), and are known to establish dorsal forebrain regional identity; (Azim et al., 2009; Bani-Yaghoub et al., 2006; Borello et al., 2014; Gaston-Massuet et al., 2016; Hagey et al., 2014; Hutton and Pevny, 2011; Johansson et al., 2013; Kikkawa et al., 2013; Kishi et al., 2012; Manuel et al., 2011; Miyoshi and Fishell, 2012; Ohtsuka et al., 2001; Ross et al., 2003; Shen and Walsh, 2005; Sur and Rubenstein, 2005; Yang et al., 2015; Zembrzycki et al., 2007).

The progenitor cell types form the triplet ${\displaystyle c_2-c_0-c_3(p\left(T=c_0|\left\{g_i^{c_0,c_2,c_3}\right\}\right)=0.96}$, Figure 4D – top right). The $c_2-c_0$ branch is mediated by 39 transition genes including LHX2, FEZF2, FOXG1, HMGA1, SP8, OTX1, SOX11, GLI3, SIX3, and ETV5 which suggest that $c_2$ is comprised of cortical progenitors (Appolloni et al., 2008; Greig et al., 2013; Kishi et al., 2012; Manuel et al., 2011; Raciti et al., 2013). The $c_0-c_3$ branch is mediated by 55 transition genes including GTF2I, HIF1A, ID1, ID3, PROX1, POU3F2, SALL1, SOX21 TCF12, TRPS1 and ZHX2, which are associated with mesencephalon and metencephalon regional development, as well as having known involvement with midbrain/hindbrain organizer identity (Buck et al., 2001; Enkhmandakh et al., 2009; Inoue et al., 2012; Jaegle et al., 2003; Kunath et al., 2002; Lavado and Oliver, 2007; Milosevic et al., 2007; Ohba et al., 2004; Uittenbogaard and Chiaramello, 2002; Yao et al., 2017). We additionally inferred the triplet ${\displaystyle c_2-c_0-c_3(p\left(T=c_0|\left\{g_i^{c_0,c_1,c_3}\right\}\right)>0.99}$, Figure 4D – bottom left), which suggests a three way split from early progenitor $c_0$ into early differentiated neuron $c_1$, and progenitors $c_2$ and $c_3$.

The continuation of the $c_2$ branch is inferred through triplet $c_0-c_2-c_4$ (${\displaystyle p\left(T=c_2|\left\{g_i^{c_0,c_2,c_4}\right\}\right)=0.97}$). The $c_2-c_4$ branch includes transition genes BCL11A, EMX2, FOXP2 and RORB, which are known to be associated with the neocortex and neuronal identity (Cánovas et al., 2015; Ebisu et al., 2016; Greig et al., 2016; Jabaudon et al., 2012; Wiegreffe et al., 2015; Woodworth et al., 2016; Zembrzycki et al., 2007).The triplet $c_0-c_3-c_5$ (${\displaystyle p\left(T=c_3|\left\{g_i^{c_0,c_3,c_5}\right\}\right)>0.99}$) meanwhile is characterized by transition genes in the $c_3-c_5$ branch including ASCL1, FOXP2, PAX3, POU3F4, ZIC1, ZIC4, HOXB2, and EN2, which have been shown to regulate fate acquisition in the midbrain/hindbrain (Agoston et al., 2012; Ang, 2006; Di Bonito et al., 2013; Elsen et al., 2008; Hegarty et al., 2013; Miller et al., 2011; Tan et al., 2014).

The $c_5$ cell cluster differentiates into two distinct clusters of post-mitotic neurons -- $c_6$ and $c_{7 }$ (${\displaystyle p\left(T=c_5|\left\{g_i^{c_5,c_6,c_7}\right\}\right)>0.99}$, Figure 4D – bottom right). The $c_5-c_6$ branch is inferred from transition genes including PAX3, CRX, SOX11, EBF2, FOXP4, ASCL1, FOXO3, and SIX3 which have strong expression in developing dopaminergic and gabaergic neurons of the midbrain (Agoston et al., 2012; Erickson et al., 2010; Pino et al., 2014; Yang et al., 2015; Yin et al., 2009; Zhang et al., 2002). Transition genes in the $c_5-c_7$ branch include ARGFX, DUXA, HES1, NFIB, PPARA, SOX2, SOX7, and SOX9, which are known to be associated with the medulla oblongata region of the hindbrain (Fawcett and Klymkowsky, 2004; Kameda et al., 2011; Kumbasar et al., 2009; Madissoon et al., 2016; Matsui et al., 2000; Stolt et al., 2003).

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, ${\displaystyle {\mathcal{𝒪}}_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.

Single-cell transcriptomic data from the in vitro neural differentiation procedure was obtained as described in Yao et al. (2017) (Supplemental information):

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 μM non-essential amino acids, 0.5 mg/mL BSA, 1X Pen-Strep, and 100 μM 2-mercaptoethanol (referred to as basal medium; all from Thermo Fisher, Waltham, MA) with 20 ng/mL FGF2 (Thermo Fisher) and 2 μM thiazovivin. Cortical induction was initiated by changing to the basal medium with 5 μM SB431542 (StemRD, Burlingame, CA), 50 nM LDN193189 (Reagents Direct, Encinitas, CA) and 1 μM 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 μM 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 μM thiazovivin. On D27, medium was changed to a 1:1 mixture of DMEM/F12 and Neurobasal medium (#21103–049) supplemented with 100 μM 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 and Yao 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 4—source data 4).

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 data 2, Figure 4—source data 5).

Calculations were performed using custom written MATLAB code (The Mathworks) on the Harvard Research Computing Odyssey cluster. Code is available at https://github.com/furchtgott/sibilant. t-SNE was done using the package provided in (Van Der Maaten, 2009). Monocole was run in R (Core Team, 2015; Trapnell et al., 2014).

The algorithm proceeds according to the following steps:

1. Find initial seed clustering configuration ${\left\{C\right\}}^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 $\left\{{\alpha }_i\right\}$ and $\left\{{\beta }_i\right\}$ given ${\left\{C\right\}}^m$:

   1. Compute ${\displaystyle p\left(g_i^{A,B,C}|T,{\alpha }_i,{\beta }_i,{\left\{C\right\}}^m\right)}$ by integrating numerically over $p\left(\mu , \sigma \right)$ (Equations 6, 9, and 10).

   2. Compute ${\displaystyle p\left(T,\left\{{\alpha }_i\right\},\left\{{\beta }_i\right\}|\left\{g_i^{A,B,C}\right\},{\left\{C\right\}}^m\right)}$ using Equations 20 and 29.

   3. Identify mostly likely topology $T$ and set of ${\displaystyle \left\{{\alpha }_i\right\}}$ and $\left\{{\beta }_i\right\}$.

3. Recluster $\left\{g\right\}$ using K-medoids in the space of all $\left\{{\alpha }_i\right\}$ and $\left\{{\beta }_i\right\}$ for the triplets with probability ${\displaystyle p\left(T|\left\{g_i^{A,B,C}\right\},{\left\{C\right\}}^m\right)}$> 0.6 of being non-null. Determine new clustering configuration ${\left\{C\right\}}^{m+1}$.

4. Repeat steps 1 and 2 until convergence of $\left\{C\right\}.$

5. Determine the most likely tree connecting cell clusters, recapitulating high-probability triplet topologies.

Each of these steps is described in the following section.