scTDA for temporal ordering of differentiating cells

We assessed the capacity of scTDA to order asynchronously differentiating cells using a controlled setting where the truth is known. First, we simulated a noisy, branched, asynchronous cellular differentiation process with two branching points (Fig. 2a and Online Methods). 700 cells were sampled at three time points and the expression levels of 500 genes were simulated. From this data, scTDA correctly reconstructed the topological structure of the differentiation tree, and identified the most stem-like state (Fig. 2b). In contrast, the algorithms Diffusion Pseudotime⁶, Wishbone⁷ and SLICER⁸ failed to correctly assign some of the branches in our hands, producing artifacts in the inferred pseudo-temporal ordering of the cells (Figs. 2c to 2e), even if the most stem-like state was provided by the user of these algorithms.

Open in a separate window

Figure 2

Comparison of algorithms for ordering cellular states

a A noisy, branched, asynchronous cellular differentiation process was simulated. The differentiation tree of the process is represented. 700 cells were sampled from this process at three time points. Using this data, we attempted to reconstruct the structure of the differentiation tree with scTDA, and compared the outcomes with reconstruction using Diffusion Pseudotime, SLICER and Wishbone, all of which rely on branch assignments for downstream statistical analysis.

b. Reconstructed differentiation trajectory using scTDA. scTDA recovered the structure of the simulated differentiation process and correctly rooted the tree using the experimental information. Nodes correspond to sets of cells sharing similar global transcriptional profiles, with the node sizes proportional to the number of cells in the node. Nodes that are connected by an edge have at least one cell in common. For reference, an inset with the latent MDS representation of the data, colored by sampling day, is also shown.

c. Reconstructed differentiation trajectory using Diffusion Pseudotime. Although the representation of the data using the first two diffusion coefficients reproduces the structure of the differentiation tree, the branches are not correctly assigned.

d. Reconstructed differentiation trajectory using SLICER. The representation constructed by SLICER using locally linear embedding was unable to capture the complete structure of the differentiation tree and branch assignments.

e. Reconstructed differentiation trajectory using Wishbone. The t-SNE representation of the data used by Wishbone reproduces the structure of the differentiation tree and correctly identifies the first branching point. However, Wishbone failed to identify the second branching point.

f, g. Pearson's correlation coefficient between the pseudo-time, inferred from the data by scTDA, Diffusion Pseudotime (D. Pt), SLICER, and Wishbone, and the actual simulated differentiation time. Cellular differentiation processes with one, two, or three branching points were considered, both in the absence (f) and the presence (g) of drop-out events.

To quantify the performance of scTDA, we extended the above simulations to processes with one, two, or three branching points (Online Methods). scTDA showed higher correlation between the inferred pseudo-time and the actual simulated differentiation time compared to the other three methods, thus reconstructing the underlying differentiation process more accurately from the data (Fig. 2f). The improvement was particularly large when drop-out events were included in the simulation (Fig. 2g), where the graphical representations produced by Diffusion Pseudotime, SLICER, and Wishbone were often unable to correctly infer the structure of the differentiation tree (Supplementary Fig. 2).

Motor neuron differentiation analysis with scTDA

We applied scTDA to longitudinal single-cell RNA-seq data from in vitro motor neuron differentiation experiments. We differentiated mESCs into motor neurons using a well-established protocol². The mESC line expresses enhanced green fluorescent protein (eGFP) under the control of the early motor neuron-specific promoter Mnx1, allowing identification of committed motor neurons. Using a modified CEL-seq approach²⁰ (Online Methods and Supplementary Note 2), we generated cDNA libraries from individual cells, spanning the transition from pluripotency to post-mitotic commitment. To assess reproducibility, we performed two separate differentiations. In the first differentiation (experiment 1) we sequenced 440 cells sampled throughout the entire course of motor neuron differentiation. Cells were concomitantly sampled in bulk at each time point. Using this data, we assessed library saturation and the optimal depth of sequencing coverage. In the second differentiation (experiment 2) we sought to increase statistical power and sequenced 2,304 individual with the depth of sequencing coverage near the saturation point (~250k reads/cell). In both experiments cells were sampled each day throughout the entire course of the differentiation, from day two to six. In total we detected 15,609 and 19,009 transcripts respectively across experiments 1 and 2. Given the improved statistical power associated with an increased number of cells, we focused our attention on the results from experiment 2, and assessed the reproducibility of the results using experiment 1.

Cell populations and surface markers of differentiation

We used scTDA to analyze the two longitudinal motor neuron differentiation single-cell RNA-seq datasets. We filtered the data for cells that passed stringent quality control tests (Online Methods, Supplementary Fig. 3), retaining, respectively, 373 and 1,964 individual cells in experiments 1 and 2. An scTDA analysis recapitulated chronological order based on mRNA expression levels alone (Figs. 3a and 3b, and Supplementary Figs. 4 and ​and5),5), and simultaneously detected detailed transcriptional relationships between individual cells. We observed only a mild dependence of library complexity on timescale of our differentiation (Supplementary Fig. 6, Online Methods), unlike other experimental settings where differences in cell cycle give rise to substantial variation of the library complexity^{21, 22}. We did not observe large batch effects (Supplementary Fig. 7). Compared with analyses of the same data sets using PCA, MDS, t-SNE, Monocle, Wishbone, Slicer and Diffusion Pseudotime, we showed that scTDA had superior ability to identify the continuous chronological structure of the differentiation process (Fig. 3c and Supplementary Figs. 8 and 9).

Open in a separate window

Figure 3

Differentiation of pluripotent mESCs into motor neurons (MN)

a scTDA recapitulates chronological order based on expression data alone. The topological representation of the expression data of 1,964 single cells, sampled from the differentiation of mESCs into MNs, is labeled by sampling time. The root node, inferred from the experimental chronological information, is indicated with a red circle.

b. The distance of each node to the root node, represented as a function of sampling time. The chronological time of a node is defined as the mean of the sampling times of the cells in the node.

c. Comparison with standard dimensional reduction algorithms. Dimensional reduction of the same expression data of 1,964 single cells, sampled from the differentiation of mESCs into MNs, using PCA, MDS, and t-SNE. The Pearson's correlation coefficient between the sampling time and the two-dimensional Euclidean distance to the root cell (defined as the one that maximizes this correlation) is indicated in each case.

d. Consistency between experiments 1 and 2. Left: Venn diagram of genes with significant gene connectivity (q < 0.05) in the topological representations of the two datasets. Both experiments are highly consistent in their calls (Fisher exact test p-value < 10^-100). The number of significant genes is larger in experiment 2, consistent with its higher statistical power (due to the larger number of cells considered). Right: correlation between the centroid (expressed in pseudo-time) of the n = 2,701 genes with significant (q <0.05) gene connectivity in both topological representations. The distribution of gene centroids is highly consistent across the two experiments.

e. Known markers of pluripotent cells, motor neuron progenitors and post-mitotic neurons are successfully ordered. The topological representation is labeled by mRNA levels of Oct3/4 (red), Olig2 (green) and VAChT (blue).

Open in a separate window

Figure 5

Differential gene expression in neurogenesis

a Differentially expressed regulators and downstream genes in the four cell populations are shown. Genes are annotated according to their role in transcriptional, cell cycle, RNA binding protein regulation and RA response.

b. Topological representations (labeled by mRNA levels) of Hoxb5, Hoxb6 and the antisense lncRNA Hoxb5os, showing concordant expression of these transcripts during the generation of motor neurons from mESCs.

c. Post-mitotic neuronal populations. Differentially expressed genes between Vsx2⁺, Gata3⁺ and En1⁺ cells that were marked as post-mitotic neurons by the scTDA analysis. Hierarchical clustering of cells leads to five groups, four of which correspond to MNs, and V1, V2a and V2b interneurons. Hierarchical clustering of genes produces four groups of genes, uniquely expressed by each of the above cell types, and a fifth group associated to GABAergic neurons. lncRNAs are marked in blue.

Analysis of the data from experiment 2 identified 7,620 genes with significant gene connectivity (q-value < 0.05, Supplementary Fig. 10 and Supplementary Table 2), compromising 74% overlap with those identified from experiment 1 (p-value < 10^-100, Fisher exact test, Fig. 3d). It is likely that the large number of significant genes is due to the transcriptional heterogeneity of the dataset, because the samples of cells represent multiple developmental stages. The centroids of the significant genes were consistent between the two datasets (Pearson correlation = 0.9, p-value < 10^-100, Fig. 3d). We validated our findings by identifying known markers of pluripotent cells, progenitor cells and mature motor neurons in the dataset over time (Fig. 3e and Supplementary Fig. 4c). We also found that cell cycle genes were downregulated (as expected) in post-mitotic neurons (Supplementary Fig. 11). We verified that the choice of parameters for building the topological representation (number of genes, bins, and overlap between bins) did not affect the result of the analysis (Supplementary Fig. 12).

Interpretation of the distribution of the centroid and the dispersion of genes (Supplementary Fig. 13) using scTDA revealed four transcriptionally distinct cellular populations arising during differentiation from mESCs to motor neurons (Fig. 4a). 488 of 7,620 genes were identified in the four expression groups (Fig. 4b, and Supplementary Table 2). Groups 1a and 1b contain genes that are only expressed in early EBs (Oct3/4⁺ cells), and correspond to pluripotent and neural precursor states. Genes in groups 2 and 3 are only expressed in the progenitor (Olig2⁺ cells) and post-mitotic ensembles (VAChT⁺ cells), respectively.

Open in a separate window

Figure 4

Cellular populations during motorneuron differentiation

a scTDA identifies four transient populations in mESC differentiation into MNs. Represented is the topological representation (colored by mRNA levels) of four groups of low-dispersion genes; pluripotent, precursor, progenitor and post-mitotic populations. In total, 488 genes were assigned to one of these four populations based on their expression profiles in the topological representation.

b. Reconstructed expression timeline for each of the four groups of low-dispersion genes.

c. Validation by detection of state-specific cell surface markers identified by scTDA. Topological representation (colored by mRNA levels) of surface proteins Pecam1, Ednrb and Slc10a4, and immunostaining of cultured EBs is shown. The scale bar in the immunostaining images denotes a length of 50 μm. Details of three regions are presented on the right. For reference, the topological representation colored by mRNA levels of the Mnx1eGFP reporter is also shown.

d. In vivo validation of the motor neuron surface marker Slc10a4. Spinal cord section of an E9.5 mouse immunostained for Slc10a4 (red). The pool of motor neurons is also marked by Mnx1eGFP expression (green). The scale bar denotes a length of 50 μm.

Ontology enrichment analyses using PANTHER²³ revealed developmental genes and genes related to DNA replication in groups 1a, 1b, and 2, whereas group 3 was enriched for genes related to axonogenesis, neuron migration, and regionalization (Supplementary Table 3), consistent with the underlying cellular differentiation process.

The analysis identified 74 transcripts in the four groups that encode proteins with an extracellular domain. These are therefore putative cell surface marker candidates that could be used to purify niche populations for further study (Supplementary Table 2). Several of these surface markers have not been reported as markers of neurogenesis. We selected three of these 74 markers that had high expression level and specific antibodies available for experimental validation. The three markers were validated using immunohistochemistry in vitro in EBs (Fig. 4c, and Supplementary Fig. 14) and in vivo in the murine embryonic spinal cord (Fig. 4d).

scTDA identifies distinct proliferative states

scTDA analysis of our motor neuron differentiation data revealed numerous loops in the neural precursor population (Fig 4a) that were not present in the mESCs or the motor neurons. These loops were also observed in the scTDA analysis of experiment 1 (Supplementary Fig. 4). Neural precursor cells proliferate rapidly upon induction with retinoic acid (RA). This is in keeping with biological expectations, given that the cell cycle will give rise to periodic structures in the expression space. To evaluate whether the loops in the topological representation of neural precursors are caused by differences in the cell cycle, we built a topological representation using only cell cycle genes (Supplementary Fig. 15a, Online Methods). The cell cycle topological representation contained substantially larger loops in the same region (Supplementary Fig. 15a), separating Stra8^up neural precursors and progenitors into proliferative and non-proliferative populations according to the expression of Mki67 and other markers of proliferation (Supplementary Figs. 15b and 15c). We used persistent homology to statistically assess the significance of these loops in the topological representation, asking whether loops of similar size could arise from noise effects (Online Methods). This analysis showed a strong statistical significance for some of the larger loops (q-value < 0.05, Supplementary Fig. 15d), consistent with a biological, rather than technical, origin for these features.

Characterization of developmental transitions

We next used scTDA to characterize the developmental transitions that occur during the differentiation of motor neurons. These results are available through an online database (Online Methods). We used scTDA to identify transcriptional programs that are associated with the transitions between the four transient cellular populations that we identified (Fig. 5a). These findings were consistent across both experiments.

First, the transition from a pluripotent cell to a neural precursor population is characterized by the transcriptional dynamics of pathways involved in RA signaling and downstream effector proteins^{24, 25} (Fig. 5a). Our analysis resolved with unprecedented resolution the timing of the transcriptional events occurring during this transition, identifying upregulation of Stra8 and downregulation of Fgf4 as some of the earliest events that mark the transition (Supplementary Fig. 16). Subsequently, there is transcriptional activation of a subset of the homeobox gene family, including Hoxa1 and Hoxb2-8, which continue to be expressed during later stages of the differentiation process, and Hoxb1, which is transiently expressed along with caudalizing transcription factors (Figs. 5a and 5b, and Supplementary Fig. 17).

A second wave of RA inducible gene activation was identified during the formation of neural progenitors. This is accompanied by transcriptional inactivation of Stra8, and activation of Crabp1 and a second set of homeobox genes, Hoxa2, Hoxa3, Hoxc5, Hoxd3 and Hoxd4 (Fig. 5a). This pattern of Hox gene activation suggests that the linear chromosomal arrangement of the Hox gene clusters does not necessitate temporal co-linearity in anterior Hox gene expression, a phenomenon that has been reported in the developing spinal cord²⁶. Both waves of homeobox gene expression were accompanied by the upregulation of several long non-coding RNAs (lncRNAs) derived from the opposite strand (Fig. 5b, Supplementary Fig. 17, and Supplementary Table 4), consistent with previously identified lncRNA-based regulation of homeobox gene clusters^27-30.

The transition to mature neurons was marked by exit from the cell cycle and post-mitotic differentiation. In keeping with expectations from previous data³¹, scTDA identified Neurog1/2 (modulators of neuronal specification, cyclin regulation, cell cycle exit) as well as Ascl1 and Sox9 among the known mediators of neuronal commitment, concomitant with a marked repression of topoisomerase 2A in the neuron population (Fig. 5a). We observed variability in eGFP expression levels and fluorescence intensities among differentiated neurons, which might indicate heterogeneity in later stages of differentiation. We performed differential expression analysis of En1⁺, Gata3⁺, Vsx2⁺ and Egfp⁺ cells with in this cellular population (Online Methods), and classified differentiated neurons into motor neurons (n = 343), and V1 (n = 19, En1⁺), V2a (n = 10, Vsx2⁺) and V2b interneurons (n = 15, Gata3⁺) (Fig. 5c), confirming the presence of cellular heterogeneity.

Interestingly, we observed the presence of several lncRNAs in the differentially expressed genes. V1 interneurons specifically express Gm12688, an intergenic lncRNA located near Foxd3, and transcribed from the opposite strand. We validated this finding using qPCR (Supplementary Fig. 18). Similarly, we identified a lncRNA in chromosome 15 (NONCODEv4 accession number NONMMUG015572) which was only expressed in V2b interneurons. Gm14204, an intergenic lncRNA located near Slc32a1 and transcribed from the opposite strand was only expressed in V1 and V2b GABAergic interneurons. These results may suggest a role for lncRNAs in neural diversification.

Among the genes identified by scTDA as being associated with developmental transitions were 20 genesencoding RNA-binding proteins (Supplementary Table 2). These include known developmental state-dependent pre-mRNA splicing factors, as well as stage-specific but uncharacterized RNA binding proteins, which may guide cellular differentiation and post-mitotic commitment. In the context of the progenitor to neuron transition, our analysis identified Nova1/2, Rbfox3, Srrm4, and the Elav-like transcripts Elavl4 and Celf3 (Fig. 5a), consistent with previous data^{32, 33}. scTDA further revealed the upregulation of Mex3b in the progenitor and post-mitotic cell populations, and constitutive expression of Ptbp1 and Ptbp2. Previously published studies have documented Srrm4 directed inclusion of neural specific exon 10 in Ptbp1³⁴. Our results may indicate a transcriptional switch in splicing factor regulation that culminates in neural specific splicing.

Single cell library generation

To obtain single cell expression profiles, modified CEL-Seq²⁰ was carried out as described in the Supplementary Note 2. Briefly, we reverse transcribed the mRNA in each cell lysate with barcoded primers, pooled the single cell samples, synthesized second strand cDNA, and then linearly amplified by in vitro transcription with T7 RNA polymerase. The amplified RNA (aRNA) was chemically fragmented and T4 RNA ligase was used to ligate Illumina 3′-RNA adapters. The aRNA was then run on an Agilent Bioanalyzer to assess proper fragmentation and then reverse transcribed to generate cDNA and subjected to PCR enrichment. The subsequent multiplexed samples were paired end sequenced (2×125 bps) on an Illumina HiSeq 2500. Bulk RNA samples were purified from 2×10⁶ cells in 1 mL of Trizol and standard Qiaquick RNA extraction protocols with a RIN 9.8 or higher. Stranded RNA-seq libraries were generated at the New York Genome Center using a TruSeq Stranded Total RNA Library Prep Kit. Stranded cDNA libraries were paired end sequenced (2×125 bps) on a HiSeq 2500, operating in high output mode, yielding 30M reads per indexed library.

Immunofluorescence

EBs were washed three times in ice cold PBS and fixed for two hours at room temperature in 4% PFA. They were then washed with ice cold PBS and either embedded into OCT and stored at -80°C for sectioning, or stored in PBS at 4°C for whole EB staining. 20um sections were cut using a Leica CM 1950 cryotome and mounted onto a glass slide. Whole and sectioned EBs were incubated in blocking solution (0.1% tween-20, 5% donkey serum, 1% BSA) for two hours at room temperature followed by primary antibodies diluted in blocking solution overnight. Primary antibodies used were: goat-αSlc9a1 (Santa Cruz sc-16097, 1:100), rabbit-αOlig2 (Millipore AB9610, 1:100), mouse-αOlig2: clone 21F1.1 (Millipore MABN50, 1:100), rabbit-αEdnrb (Novus NLS54, 1:200), guinea pig-αVAChT (gift from Dr. Neil Schneider's lab, 1:400), rabbit-αSlc10a4 (Novus Biologicals NBP1-81134, 1:100), mouse-αPou5f1 (BD Biosciences 611202, 1:100), goat-αPECAM1 (Santa Cruz sc-1506, 1:100), sheep-αCD44 (R&D Systems AF6127, 1:250), goat-αFoxc1 (Santa Cruz sc-21396, 1:50), and goat-αSstr2 (Santa Cruz sc-11606, 1:50). Alexa Fluor-conjugated secondary antibodies from Life Technologies were used at a 1:1000 dilution for two hours at room temperature. CD44, Foxc1 and Sstr2 were respectively conjugated to Cy-5, Alexa-405, and Cy-3 fluorophores using the DyLight system from Abcam (ab201798, ab188287, ab188288). Coverslips were mounted with Vectashield and EBs were imaged on an Olympus Fluoview FV1000 microscope using Olympus Fluoview v4.1.

Processing of single cell RNA-seq data

Paired-end 125bp reads were de-multiplexed, trimmed, and mapped to the UCSC mouse reference (mm10) using Tophat³⁸. Gene expression was quantified using transcript read counts as derived from HTSeq³⁹. Read counts were normalized as,

where r_i denotes the unambiguous read count for transcript i, and r denotes the total number of reads that are mapped to transcripts in the cell. A strategy based on spike-in read counts, as described by Stegle et al.⁴⁰, was implemented to filter out cells with a low content of mapped RNA and/or low sequencing depth. Specifically, the ratio ₁ between average spike-in reads in the cell and average spike-in reads in the library (for spike-ins with an average of more than five reads in the library) was used to discard cells with low sequencing coverage. Similarly, the ratio ₂ between the total number of spike-in reads in the cell and the total number of mapped reads was used to discard cells with a low number of mapped reads, relative to the sequencing coverage. Those cells showed very low expression of house-keeping genes, possibly representing cells under stressed conditions and/or with large amounts of degraded RNA. Based on the distribution of ₁ and ₂ in each of the two experiments (Supplementary Fig. 2a), cells with ₁ > 5 · 10⁻⁴ and 0.7 > ₂ > 0.01 in the experiment 1, and cells with 4.0 > ₁ > 0.05 and 1.0 > ₂ in experiment 2, were kept for subsequent analysis. The distribution of filtered out cells across different libraries was uniform in both experiments (Supplementary Fig. 2b). To reduce the noise observed near detection threshold in experiment 1, read counts with r_i < 5 were set to zero. Additionally, one of the libraries (RPI36) was discarded from subsequent analysis, because it presented a large batch effect. To assess the dependence of the library complexity on the differentiation time (Supplementary Fig. 3), we computed at each time point the distribution of the geometric library size, defined as the sum of log expression values over all genes in a cell²².

Topological representation

The algorithm Mapper¹⁹ (Supplementary Note 1) was used to build topological representations of the RNA-seq data, through the implementation provided by Ayasdi Inc. Several open-source implementations of Mapper are available (https://github.com/MLWave/kepler-mapper, http://danifold.net/mapper/, https://github.com/RabadanLab/sakmapper, https://github.com/paultpearson/TDAmapper). In brief, the processed RNA-seq data was endowed with a dissimilarity matrix by taking pairwise correlation distance (1 - Pearson correlation). To minimize the effect of drop-out events present in single-cell data, we only considered the 5,000 genes (for experiment 1) and the 4,600 genes (for experiment 2) with highest variance across each dataset. These are highly expressed transcripts for which the probability of not being captured by the RNA amplification (drop-out events) is small^41-44. The space was reduced to R² using MDS, as displayed in Fig. 3c. A covering of R² consisting of 26×26 (62×62) rectangular patches was considered for experiment 1 (2), respectively. The size of the patches was chosen such that the number of cells in each row or column of patches was the same, hence avoiding sampling density biases. The overlap between patches was 66% (on average). Single linkage clustering was performed in each of the pre-images of the patches using the algorithm described in Singh et al.¹⁹. A network was constructed in which each vertex corresponds to a cluster, and edges correspond to non-vanishing intersections between clusters. We checked for the absence of batch effects in the topological representation by verifying that batches from the same timepoint had a substantial overlap in the representation (Supplementary Figs. 3d and 7) and the stability against different choices for the threshold for the number of genes used to compute the distance matrix and for the covering of R² (Supplementary Fig. 12).

We verified that our analysis was reproducible using a freely-available open-source implementation of Mapper. To that end, we repeated the analysis of experiment 1 using a topological representation built with Sakmapper (https://github.com/RabadanLab/sakmapper) (Supplementary Fig. 20). Comparison with the results of the scTDA analysis presented in the main text shows a clear consistency between the two analyses. In particular, 78% of the genes with significant gene connectivity (q-value < 0.05) in the representation generated with Sakmapper also displayed significant gene connectivity in the Ayasdi representation (Fisher exact test p-value < 10^-100). The centroid and dispersion of those genes were also consistent between the two topological representations, with Pearson's correlation coefficient r = 0.98 (p-value < 10^-100) for the centroid (Supplementary Fig. 20b), and r = 0.94 (p-value < 10^-100) for the dispersion (Supplementary Fig. 20c) of significant genes.

Gene connectivity, centroid and dispersion within the topological representation

A notion of gene connectivity in the topological representation was introduced, defined as

where e_i,α represents the average expression of gene i in node α of the topological representation, normalized as described in the paragraph “Processing of single cell RNA-seq data”; Γ denotes the set of nodes of the topological representation; A_αβ is its adjacency matrix; and N the total number of nodes. With this normalization, s_i takes values between 0 and 1.

The gene connectivity score depends on the distribution of expression values of the specific gene (Supplementary Fig. 10), and therefore genes cannot be ranked accordingly to their gene connectivity score in a meaningful way. To assess the magnitude of the connectivity score relative to genes with the same expression profile, and rank genes accordingly, we introduced a non-parametric statistical test. We tested for the null hypothesis of a randomly expressed gene with the same distribution of expression values having a higher gene connectivity score. To that end, a null-distribution was built for each gene i using a permutation test. Cell labels were randomly permuted 5,000 times for each gene, computing s_i after each permutation. A p-value was estimated by counting the fraction of permutations that led to a larger value of s_i than the original one. Gene connectivity and its statistical significance were computed for each gene expressed in at least three cells. The resulting p-values were adjusted for multiple testing by using Benjamini-Hochberg procedure for controlling the false discovery rate.

To establish a pseudo-temporal ordering within the topological representation, a notion of root node was introduced. The latter was defined as the node that maximizes the function,

where corr(x,y) denotes Pearson's correlation coefficient between x and y; d_α is the graph distance function to node α, that assigns to each node of the topological representation a value corresponding to the number of edges that are crossed in the shortest path from that node to node α; and t is the chronological sampling time function, that assigns to each node of the topological representation the average sampling time (expressed in days) of the cells contained in the node.

Least-squares linear regression was performed to determine the best fit for the coefficients a₀ and a₁ in the relation

where d_root is the graph distance function to the root node, determined in the previous paragraph. These coefficients were used to define the centroid and dispersion of each gene in the topological representation, expressed in days, and given respectively by

and

Such normalization, i.e. using coefficients a₀ and a₁ to express the centroid and dispersion in units of pseudo-time (days), allows comparing the connectivity and dispersion of a gene across different topological representations or studies.

Significance of topological features

We computed the first persistent homology group^{45, 46} using the graph distance of the topological representation. Given the pairwise distances of a set of points sampled from a space, persistent homology enables the quantification of topological features (connected components, loops, cavities, etc., preserved under continuous deformations of the space) compatible with the data at each scale. The first homology group, in particular, classifies loops of the space (Supplementary Note 1). We used persistent homology death times as a proxy of the size of the loops, and evaluated their statistical significance using a permutation test. To that end, we randomly permuted the labels of the genes 500 times, for each cell independently. For each permutation we built a topological representation using the same parameters as in the original representation and computed the first persistent homology group. A p-value for each of the loops was estimated from the distribution of the number of loops as a function of their death time. The resulting p-values were adjusted for multiple testing by using Benjamini-Hochberg procedure for controlling the false discovery rate.

Comparison with other methods for analyzing longitudinal single-cell RNA-seq data

We dimensionally reduced the processed single cell RNA-seq data of experiment 2 using MDS, ICA and tSNE. For MDS and ICA we used the same set of highly expressed, highly variant genes that we used for building the topological representation. For tSNE, we used the top 10 principal components of the gene space. In each representation we determined the cell that maximized the Pearson correlation coefficient between the two-dimensional Euclidean distance to the cell and chronological sampling time, corresponding to the least differentiated cellular state.

Additionally, we compared scTDA to the single cell software Monocle¹⁰, based on ICA and minimum-spanning trees, Diffusion Pseudotime⁶ and Wishbone⁷, based on diffusion coefficients, and SLICER⁸, based on locally linear embedding. We followed all recommendations in the documentation of these algorithms. In our tests, Monocle failed in running over the complete dataset from experiment 2, consisting of 2,304 cells, and only a partial set of 834 cells, sampled from all time points was analyzed.

Simulated data

Noisy, branched asynchronous cellular differentiation processes were simulated, from which 700 cell were sampled at three time points. To that end, we used the following strategy:

1. - We simulated a noisy branched tree-like structure with three parameters (t, u, and v), and performed a non-linear transformation into four-dimensional embedding space (spanned by the variables g⁽¹⁾, g⁽²⁾, g⁽³⁾, and g⁽⁴⁾). This space provides the structure for four different groups of correlated genes.

2. - We randomly sampled 700 points from this space, corresponding to 700 cells, and assigned a sampling day based on a multinomial distribution with probabilities given by a logistic function of t, to simulate asynchrony.

3. - We simulated the expression of 300 genes driven by the variables g⁽¹⁾, g⁽²⁾, g⁽³⁾, and g⁽⁴⁾, in addition to 200 genes with non-correlated expression, sampled from normal distributions.

4. - We simulated the effect of drop-out events using the standard logistic dependence of the drop-out probability as a function of expression.

In what follows, we provide a detailed description of each of these steps.

First, we simulated four groups of correlated genes. These were defined by the equations,

where $N_i^{\left(k\right)}$, are normally-distributed random variables with mean μ = 150 and standard deviation σ = 8. The index i = 1, …, 700 runs across the sampled cells, and u_i, v_i, and t_i are randomly sampled from

where only the first 2 + r equations are considered in simulated differentiation processes with r lineage branching points. The variable t_i represents the differentiation pseudo-time of cell i at the time of sampling. To simulate asynchrony, each sampled cell was assigned a sampling day accordingly to

where p_i [0,1] is a random number uniformly distributed, and τ is the logistic function

We simulated 75 genes in each of the four groups of genes, with expression values given by

where ${N^{\prime }}_l^{\left(k\right)}$ and ${N^{″}}_i^{\left(k\right)}$ are normally-distributed random variables with mean μ = 1 and standard deviation σ = 0.2. In addition, we simulated 200 extra genes with non-correlated expression

with ${N^{\prime }}_l^{\left(0\right)}$ and ${N^{″}}_l^{\left(0\right)}$ normally-distributed random variables with mean μ = 200 and 1, and standard deviation σ = 50 and 0.2, respectively. Hence, a total of 500 genes were simulated in each of the 700 cells. To model the effect of drop-out events, we randomly set to zero the expression of some of the genes in some of the cells, with probability

where m is the original expression value of the gene in the cell.

Apart from scTDA, we ran the algorithms Diffusion Pseudotime⁶, Wishbone⁷, and SLICER⁸ in the simulated datasets. We followed all the recommendations in the documentation of these algorithms. We assessed the performance of each algorithm by computing the Pearson's correlation coefficient between the cell pseudotime inferred by each of these algorithms and the simulated variable t.

Gene ontology annotation

Gene ontologies were obtained from EMBL-EBI QuickGO⁴⁷. Specifically, categories GO:0006355 “Regulation of transcription, DNA-templated”, GO:0008380 “RNA splicing”, GO:0044822 “Poly(A) RNA binding”, GO:0051726 “Regulation of cell cycle”, and GO:0007049 “Cell cycle” were used to annotate genes. Expression of genes associated to DNA replication was analyzed by considering the 99 genes in the category GO:0006260 “DNA replication” expressed in less than 1,400 cells in experiment 2. Genes coding for proteins in the cellular surface were identified by looking in UniProtKB database⁴⁸ for proteins annotated with an extracellular topological domain. Gene ontology enrichment analysis was performed using PANTHER classification system²³.

Transient cellular populations

Low-dispersion genes (k_i < 1.7 days and k_i < 2.25 days respectively in experiments 2 and 1) with significant gene connectivity (q < 0.05) in the topological representation were clustered according to their centroid using k-means clustering (Supplementary Fig. 13b and 13c). The optimal number of clusters according to Davies-Bouldin index was four in experiment 2 (three in experiment 1) (Supplementary Fig. 13a), as it was also evidenced from visual inspection of the centroid distribution for low-dispersion genes. A state r = 1, …, 4 was assigned to each node of the topological representation based on the average expression of each cluster of low-dispersion genes in the node (Supplementary Fig. 13d). Genes with significant gene connectivity according to the permutation test described in the paragraph “Gene connectivity, centroid, and dispersion within the topological representation” were assigned to each of the four populations based on the number of cells expressing the gene in each state r. Only genes expressed in at least 80 cells and at most 1,500 cells were considered.

Analysis of long non-coding RNAs

The coordinates of intergenic and antisense lncRNAs were downloaded from NONCODEv4⁴⁹, and read counts were obtained using HTSeq (http://www-huber.embl.de/users/anders/HTSeq). The connectivity and statistical significance of each long non-coding gene in the topological representation was computed using scTDA. Only lncRNAs that were significant (q < 0.05) in both experiments 1 and 2, and that were supported by at least 50 reads in the longitudinal stranded RNA-seq data were kept. Curation was performed to remove lncRNAs whose 3′-end overlapped the 3′-end of another gene, gene assignment therefore being ambiguous, or that corresponded to possible miss-annotations of the 3′ UTR of a nearby gene.

Characterization of interneuron populations

Differential expression analysis between En1⁺, Gata3⁺, Vsx2⁺ and Egfp⁺ cells in nodes characterized as post-mitotic (Supplementary Fig. 13d, state 3) was performed using the software SCDE⁴² with default parameters.

scTDA software

The algorithms described in this work were implemented and documented in an object oriented python library for topological data analysis of high-throughput longitudinal single-cell RNA-seq data, called scTDA. scTDA is provided as Supplementary Code and is publicly available at https://github.com/RabadanLab/scTDA.

We thank Tom Jessell, Nicole Francis, and Hemali Phatnani for critical reading of the manuscript. A.H.R. and T.M. thank the New York Genome Center and D. Goldstein for sequencing support, Susan Morton for providing Engrailed antibody and P. Sims for experimental discussions. P.G.C. and R.R. thank A. J. Levine, F. Abate, I. Filip, S. Zairis, U. Rubin, and P. van Nieuwenhuizen for useful comments and discussions, and O. T. Elliott for technical support with the online database. The work of P.G.C and R.R. is supported by the NIH grants U54-CA193313-01 and R01GM117591. The work of A.H.R., E.K.K., T.J.R. and T.M. is supported by ALS Therapy Alliance grant ATA-2013-F-056 and NIH grant NS088992.