Defining cell types and states with single-cell genomics

Defining cell types and states requires single-cell assays

Single-cell measurements help overcome several key obstacles that have frustrated efforts to precisely define cellular states and catalog them in development and disease. Except in rare instances where cells can be precisely synchronized, bulk measurements destroy crucial information by averaging signals from individual cells together. One example of this is described by Simpson's Paradox (Simpson 1951), which is well known in statistics but rarely discussed in cell biology. Consider an experiment aiming to assess whether two genes show correlated expression levels across a population of cells. Finding such pairs of genes is often a major goal of single-cell genomics experiments, as some cell fate decisions are governed by a single pair of mutually exclusive regulators. For example, the myeloid/erythroid fate decision in common myeloid progenitor cells is made via the mutually antagonistic transcription factors PU.1 (encoded by SPI1) and GATA1 (Rekhtman et al. 1999). In the hypothetical experiment, one might observe the pattern on the left side of Figure 1A and conclude that the genes in question are expressed mutually exclusively. However, if the population consisted of a mixture of two separate groups of cells, this conclusion might be incorrect. Grouping the cells properly by type, and then performing the analysis could reveal the genes are in fact positively correlated, not negatively correlated. That is, failing to properly compartmentalize the data by cell type leads to a qualitatively incorrect interpretation. The misleading effects of Simpson's Paradox are likely to be pervasive in modern experiments using bulk assays.

Another crucial reason that single-cell measurements are necessary to define cell states is that bulk measurements confound changes due to gene regulation with those due to shifts in cell type composition. Consider an experiment aimed at studying the effect of a drug on a tissue composed of two cell types. (Fig. 1B) Suppose a certain gene's expression is measured via bulk-cell analysis before and after treatment with the drug. Upon measuring a major increase in the gene's expression, one might surmise that the drug causes an up-regulation of the gene in the relevant cells. However, it might instead be that the drug stimulates cell division in one of the two subpopulations. Bulk assays cannot discriminate between these two cases.

Time series studies of gene expression, which are foundational for examining dynamic processes in development and disease progression, also suffer from averaging in bulk studies. For example, differentiating cells might transition through a sequence of intermediate states on the way to becoming fully mature. However, they typically do not do so in a synchronized manner, so sampling a population of them at any particular moment in time will yield a mixture of cells from different stages, as illustrated in Figure 1C. Tracking the time-average expression of a gene of interest might yield a qualitatively misleading picture of that gene's developmental regulation. If we could reorder the cells on the horizontal axis, we would see a faithful representation of the gene's expression dynamics during development.

Bulk measurements are thus fundamentally constrained by averaging. Accurately defining the cell types and states in our bodies and explaining how they arise in development and disease demands single-cell measurements.

Technological advances in cellular state measurement

Next-generation sequencing has proven to be a remarkably sensitive means of monitoring gene expression, epigenetic configuration, nuclear structure, and other aspects of cellular state. However, most assays require a minimum level of input material that exceeds that of a single cell. Over the past several years, a number of sequencing-based assays have been optimized to work at the level of individual cells. These improvements have primarily stemmed from breakthroughs in amplification techniques, modifications to reverse transcriptase that improve processivity and enable controlled template switching, and the development of instruments that physically capture and isolate individual cells.

Most single-cell genomics assays have been adapted from similar techniques developed for analyzing bulk-cell populations (Table 1). The most widely used of these is RNA-seq, which measures global gene expression by reverse transcribing RNA into cDNA and sequencing it (Cloonan et al. 2008; Mortazavi et al. 2008; Nagalakshmi et al. 2008). Counting the reads that originate from each gene yields a measure of its expression. Single-cell versions of the RNA-seq protocol isolate individual cells in microfluidic capillaries (Wu et al. 2013), by serial dilution, or via flow sorting into multiwell plates. Each cell can thereafter be lysed to recover its RNA. Generating cDNA libraries from picograms of RNA has required major technical innovations and remains an area of considerable development effort (Tang et al. 2009; Islam et al. 2011; Hashimshony et al. 2012; Ramsköld et al. 2012; Picelli et al. 2014). A recent review provides a thorough discussion of the molecular biology underlying recent protocols (Saliba et al. 2014).

Analyzing single-cell RNA-seq data is substantially more difficult than bulk experiments. The size and scale of the data puts considerable strain on most algorithms: Where bulk experiments typically have a few dozen samples at most, typical single-cell studies capture hundreds (Shalek et al. 2013; Trapnell et al. 2014; Treutlein et al. 2014) or even thousands (Shalek et al. 2014; Fan et al. 2015; Klein et al. 2015; Macosko et al. 2015) of cells. Furthermore, single-cell expression measurements are often highly variable, and separating technical from biological variability is essential. New normalization techniques based on exogenous “spike-in” standards or molecular barcoding (Fu et al. 2014) may be needed to attribute changes in expression across cell populations to genuine biology (Brennecke et al. 2013; Grün et al. 2014). For a detailed discussion of these issues, see Stegle et al. (2015); Kolodziejczyk et al. (2015). Finally, single-cell RNA-seq enables a biologist to ask and answer questions that simply cannot be approached with bulk data and for which no tool currently exists. New algorithms and software tools may be required, both to obtain accurate results and also to exploit the unique possibilities of these measurements.

Integrating imaging with single-cell genomics will help connect cell form, function, and communication with gene regulation. Remarkably, individual cells maintain similar concentration of individual mRNA species despite wide variation in cell volume, suggesting that global transcriptional regulation is tightly connected to the physical properties of the cell (Padovan-Merhar et al. 2015). Mats Nilsson and colleagues used rolling circle amplification of padlock probes to resolve subcellular localization of multiple highly similar RNA species in individual cells (Larsson et al. 2010). More recently, Lee et al. (2014) devised FISSEQ, which obtains mRNA abundances and subcellular localization for thousands of genes by performing cDNA sequencing-by-synthesis reactions within single cells rather than on a sequencer's glass flowcell. SeqFISH and MERFISH are combinatorial fluorescence in situ hybridization strategies that use successive rounds of labeling, imaging, and photobleaching to recover mRNA abundances and subcellular localization simultaneously (Lubeck and Cai 2012; Chen et al. 2015). These techniques could in principle illuminate gene networks that mediate cell–cell communication in solid tissues. Several groups have reconstructed spatial position of individual cells directly from their transcriptomes via statistical inference in developing mouse and zebrafish embryos (Durruthy-Durruthy et al. 2014; Achim et al. 2015; Satija et al. 2015). The intuition behind such algorithms is that cells from the same region of a tissue are receiving similar paracrine signals and reside at similar positions within morphogen gradients, which should induce a common signature in these cells. An algorithm that can identify those signatures should be able to orient all cells with respect to one another. Single-cell analysis is already helping to map novel cell–cell interactions under controlled, in vitro settings (Shalek et al. 2014). Single-cell genomic analysis of such interactions in vivo will provide a radical advance in our understanding of how cells collaborate in solid tissues.

RNA-seq is only one of several aspects of cell state to be measured with single-cell genomics. Peter Fraser and colleagues adapted Hi-C, a technique for measuring chromatin conformation (Lieberman-Aiden et al. 2009) to work in single cells (Nagano et al. 2013). Hi-C measures physical proximity between each pair of sites in the genome in three dimensions, producing a “contact map” that can be used to identify looping interactions between regulatory elements and gene loci. Single-cell Hi-C might be used to investigate a key question in development, i.e., does the highly structured packing of DNA into the nucleus, which is disrupted as part of every cell division, contribute to cell fate decisions? Do stable cellular states have corresponding “signature” chromatin conformations? Single-cell Hi-C has so far only been applied to a small number of cells, but with enough observations, it might be possible to answers such questions.

If RNA-seq measures the “output” of each gene locus, epigenetic sequencing assays such as bisulfite sequencing and DNase I hypersensitivity sequencing track a locus’s “input” signals. Whole-genome bisulfite sequencing originally introduced to assay DNA methylation (Lister et al. 2008) has been successfully adapted to work in single cells (Guo et al. 2013) through “reduced representation” mapping (Meissner et al. 2008), which assays DNA methylation at CpG islands. This snapshot of the methylome constitutes an epigenetic signature that distinguishes different gene regulatory states. Histone modifications, such as H3K4 trimethylation, which marks active gene promoters, can be assayed genome-wide by performing chromatin immunoprecipitation and then sequencing the recovered DNA (ChIP-seq). ChIP-seq for chromatin marks or transcription factors can be used to create highly informative descriptions of the gene regulatory state of a cell population (Johnson et al. 2007; Mikkelsen et al. 2007). Although ChIP-seq has not yet been adapted to work in single cells, it is being optimized for smaller and smaller input (Lara-Astiaso et al. 2014).

Just as each cell type displays characteristic gene expression programming, cell types are identified by distinct sets of genomic regulatory sequences that are accessible to DNA binding profiles (Stergachis et al. 2013). There are now several genome-wide assays for accessibility (Crawford et al. 2006; Giresi et al. 2007; Hesselberth et al. 2009; Buenrostro et al. 2013). We recently adapted ATAC-seq, which uses an engineered transposase to construct a sequencing library from accessible DNA elements (Buenrostro et al. 2013), for use in single cells (Buenrostro et al. 2015; Cusanovich et al. 2015). Tn5 not only cleaves DNA at nucleosome-free regions, it also inserts DNA sequencing library adapters in situ. Thus, treating chromatin with Tn5 produces a library of fragments that mostly come from nucleosome free regions, albeit with lower genomic resolution than DNase I hypersensitivity sequencing (DNase-seq). Remarkably, the Tn5 reaction can take place within intact nuclei, allowing them to be barcoded for single-cell sequencing without ever physically isolating individual nuclei. Physical isolation is currently a requirement of most single-cell genomics protocols and generally results in cost and effort that scales linearly with the number of cells. In contrast, cost and effort for the “cellular indexing” strategy introduced by Cusanovich et al. (2015) scale sublinearly. This scalability enables the profiling of hundreds, thousands, or in principle, tens of thousands of cells in a single experiment.

Assaying gene expression, epigenetic state, nuclear structure, and other aspects of individual cells has enormous advantages over bulk assays, but it also comes with risks and challenges that are poorly understood. One basic difficulty stems from disaggregating a solid tissue into a suspension of single cells. For example, tissue samples are typically treated with enzymes such as collagenase to degrade proteins that hold the cells together. Experiments that use single-cell genomics to survey cell types and frequencies require that the disaggregation be unbiased. A procedure that leads to lysis of one cell type in the tissue but not another would dramatically skew the results of a single-cell genomics analysis. Furthermore, the very act of dissociation could substantially perturb the state of the cells. For example, cells that rely on integrin signaling or communication via elements of the extracellular matrix might be greatly disturbed between disaggregation and measurement. Careful comparisons with bulk assays will be needed to determine the extent of bias introduce by these and other problems that arise from the need to “singularize” tissues.

Ultimately, single-cell measurements of transcriptome and epigenome state will help construct a global view of gene regulation in development and disease. What if we could profile both the transcriptome and the epigenome in an individual cell? Such measurements might unlock the “histone code,” revealing how histone modification regulates transcription (Strahl and Allis 2000). However, integrating several types of single-cell measurement poses an enormous experimental and analytic challenge. The amount of DNA and RNA in a single cell is minute: How much might one lose by splitting it among multiple assays? Nevertheless, two recent studies demonstrated how single cells can be genotyped and profiled for global gene expression simultaneously (Dey et al. 2015; Macaulay et al. 2015). Cocapture of different layers of gene regulatory information in individual cells will power statistical models of how upstream regulatory activity corresponds to transcriptional output.

A dynamical systems view of the cell

In 1957, C.H. Waddington proposed a powerful metaphor for how differentiated cell types emerge from a single, totipotent cell (Waddington 1957). The “epigenetic landscape” envisions a ball rolling down a hill, encountering ridges and furrows in the hill as it goes. These features restrict its progress, ultimately forcing it to wind up at or near one of several locations at the base of the hill. Waddington used this metaphor to provide an intuition for how the network of biochemical and gene regulatory interactions in our cells endows them with both robustness and plasticity. Waddington was inspired by dynamical systems, which are mathematical frameworks for predicting the qualitative behavior of complex phenomena, such as the orbits of planets, over time. The metaphor has been refined and formalized over time to explicitly describe a cell's state at a given time as a point in “gene regulation space.” As the cell regulates gene expression during development, it moves through the space along a trajectory. The shape of the trajectory is a function of (1) the cell's starting position in the space, and (2) the differential equations that dictate changes in the cell's gene expression levels that will occur at the next moment in time, given where it currently is in the space. The central question posed by Waddington is: Where will the cell end up?

The landscape analogy is appealing, but we are a long way from using dynamical systems to describe cell fate transitions in most developmental contexts. Nevertheless, some important steps have been taken on the road to constructing useful models of cell state transitions. Stuart Kauffmann theorized that not only can cell fate transitions be described by a dynamical system, but that what we think of as distinct cell types correspond to attractors of that system (Kauffman 1993). An attractor is a location in the state space of a dynamical system that will move toward and to which it will return if perturbed. That is, trajectories flow toward attractors. Later work by Sui Huang showed that a simple regulatory network with just two genes that inhibit one another but promote their own expression is sufficient to generate a hypothetical landscape with three stable cell fates (Huang et al. 2007). Two of these states are highly stable, but the third is only weakly stable, and a small perturbation would induce cells in this state to move to one of the other two (Fig. 2A). This situation is reminiscent of a stem cell, which is poised to choose one of several differentiated fates, but will do so only upon some sort of stimulus. Furthermore, modifying the network by weakening the tendency of the genes to promote their own expression makes the “stem-like” state less and less stable, until finally it disappears. Mathematical models describing lateral inhibition mediated by Notch-Delta interactions and the resulting changes in cell state have been experimentally validated and improved to include more members of the pathway (Collier et al. 1996; Sprinzak et al. 2010; Boareto et al. 2015). James Ferrell used nonlinear dynamics to model mitotic oscillators and examined bistability in Xenopus oocyte maturation (Ferrell and Machleder 1998). Thus, although dynamical systems models have been restricted to just a few genes, they have captured remarkably rich behaviors commonly observed in development.

Single-cell clustering and trajectory analysis

Single-cell genomics offers a means of precisely quantifying the state of individual cells and thus may enable the construction of explicit, genome-scale dynamical cellular models. Early single-cell transcriptomic studies lend support to the idea that cells are occupants of a vast, complex landscape of possible states and raise doubts that cell types are precisely defined, discrete entities (Buganim et al. 2012; Moignard et al. 2013, 2015; Xue et al. 2013; Durruthy-Durruthy et al. 2014; Kumar et al. 2014; Pollen et al. 2014; Shalek et al. 2014; Trapnell et al. 2014; Treutlein et al. 2014; Buettner et al. 2015). Even within a single cell type, there is considerable variation in global gene expression, consistent with the idea that most cells reside within wells of attraction. Time series experiments of differentiation have observed cells transitioning between a starting state and one or more end states, with many cells distributed along a trajectory between them (Bendall et al. 2014; Trapnell et al. 2014; Buettner et al. 2015; Moignard et al. 2015).

Unsupervised clustering might be an effective way to identify attractors and thus define cell types and stable states using single-cell genomics data. For example, consider an RNA-seq experiment that has captured a mixture of cells from a solid tumor. The tumor might contain a mixture of three types of cells: invasive tumor cells, noninvasive cells, and stromal cells. We would like to characterize the expression profile of each and assign the cells to the proper type. Each cell captured in the experiment can be represented as a point in a high-dimensional expression space; that is, if there are 30,000 genes in the genome, each cell is a point in a 30,000-dimensional space. Clustering the cells amounts to measuring the distances between all pairs of points, then grouping them into neighborhoods based on mutual proximity. Although there are many algorithms for solving this task, they all face a sobering mathematical reality: the “curse of dimensionality.” Distances between points become more and more similar as the dimension of the space they reside within increases. Specifically, the distance between a point and its most distant neighbor approaches the distance between the point and its nearest neighbor (Beyer et al. 1999). If all the cells are equidistant from one another, how are we to cluster them?

Fortunately, in most biological systems, we do not have to shoulder the full burden posed by 30,000-dimensional measurements. Most genes can be grouped together into “modules” with expression levels that are highly correlated across cells in the experiment. Although a module might contain 1000 genes, if they are tightly correlated, we only need one representative gene to describe them all. The cells in a single-cell RNA-seq experiment can often be represented in a space with far fewer than 30,000 dimensions: We only need one dimension for each module of genes. As for unsupervised clustering, there are many algorithms that aim to reduce the dimensionality of data. Principal components analysis (PCA) is by far the most widely used, and it has already proven effective in a number of single-cell genomics studies. By “projecting” the cells down to the first two principal components, we are describing each cell as a point in two dimensions instead of 30,000. For example, several groups have studied developing embryos (Hashimshony et al. 2012; Grindberg et al. 2013; Xue et al. 2013) with single-cell transcriptomics by projecting them onto the first two principal components. In PCA space, cells from two-cell embryos are close to one another and closer to those from four-cell embryos than eight-cell embryos, with later stages still further away. Are two dimensions sufficient to capture all of the subtleties and nuances contained in a real experiment? Probably not, but reducing dimensionality might enable grouping the cells by type, which might be required for other downstream analyses.

If clustering single cells reveals the attractors of the cell state space, can it also be used to find transition paths between stable states? Can we reconstruct cellular trajectories that correspond, for example, to cell differentiation, progression through the cell cycle, or other dynamic processes? Recently, several groups have developed algorithms for not only clustering cells by type, but also organizing them by temporal or developmental stage. These algorithms have aims similar to prior strategies (Magwene et al. 2003; Qiu et al. 2011) for staging embryos, tumors, and other settings, where progression is not a simple function of time. The Monocle algorithm introduced the notion of pseudotime, a quantitative measure of biological progression through a process such as cell differentiation. Monocle first reduces the dimensionality of the expression data, then reconstructs a trajectory along which the cells are presumed to travel, and finally projects each cell onto this trajectory at the proper position. Each cell's pseudotime value is measured as the distance along the trajectory from its position back to the beginning (Fig. 2B). In order to describe complex differentiation processes in which cells make fate decisions, Monocle allows these trajectories to have a branched structure with multiple possible outcomes or “lineages” (Trapnell et al. 2014). Monocle is designed to work with single-cell RNA-seq data, but a similar approach for proteomic data was recently described by Dana Pe'er and colleagues, whose algorithm, Wanderlust, organizes cells measured with high dimensional cytometry (CyTOF) using an entirely different algorithm than Monocle, and achieves similar ends: Bendall et al. (2014) captured millions of differentiating B cells and correctly ordered them by a proteomic measure of pseudotime. More recently, Moignard et al. (2015) used single-cell qPCR to capture the bifurcation between blood-forming and endothelial cells with yet a third algorithm based on diffusion maps. These algorithms all assume that cells are sampled at different points along a complex biological trajectory. Regulation of the cell cycle is intimately coupled to cell differentiation, which could confound trajectory analysis. Should the signal driven by the cell cycle be removed from the data to bring differentiation-dependent changes into focus? Should it be explicitly modeled, even if that means sacrificing power in downstream analyses? Buettner et al. (2015) introduced an algorithm that captures the signal due to the cell cycle (or in principle, other “confounding” phenomena) that is compatible with trajectory analysis. Although these techniques are still very new and a great deal of work remains to be done, it is clear that single-cell experiments can illuminate trajectories in a wide variety of biological settings.

Finding genes that govern cell state

Ultimately, most single-cell studies aim to identify and characterize the genes that determine where a cell is in the state space and drive a cell to transition from one state to another. Doing so requires first distinguishing the genes that are differentially expressed between states or during state transitions from genes that are constitutively expressed throughout the state space. Differentially regulated genes must then be further analyzed to identify those that are responsible for establishing cell state or drive transitions. Solving these two problems may sound straightforward. However, both are surprisingly difficult (the second far more than the first) and have been the subject of countless studies.

Differential expression analysis refers broadly to the task of identifying those genes with expression levels that depend on some variable, like cell type or state, time since perturbation, or genetic background. A biologist that wishes to identify such genes in an experiment collects RNA samples from two or more conditions and compares them with a program designed to find those that differ more than some critical threshold. Exactly what that threshold needs to be depends on how much measurement variability exists in the assay as well as how much uninteresting biological variability exists within each condition in the contrast. To control for uninteresting changes in expression, each condition can be replicated, either by measuring a sample multiple times (to assess technical variability) or repeating the experiment several times (to assess biological variability). The algorithm used to analyze the expression data must then learn a statistical model that accounts for this variability and use it to set the critical threshold for changes in gene expression. This problem has been extensively studied and is essentially solved for expression microarrays, conventional RNA-seq, and a number of other bulk assays.

Single-cell genomics demand new algorithms and software for identifying genes that are differentially expressed. Current single-cell protocols destroy a cell in order to generate a library from it, so a given cell can only be measured once. Although in principle, the lysate from a single cell could be split among several technical replicates, this capability has not yet been demonstrated. Thus, in single-cell genomics experiments, there are no “replicates” per se, which increases the risk that studies will be confounded by batch effects. As discussed above, controlling for library quality and accounting for technical variability is a major focus of technological developments. Analysis software will need to explicitly accommodate protocol improvements in order to distinguish biological variation from technical noise across cells. Single-cell experiments also typically have far more samples than bulk experiments. Due to cost and labor concerns, bulk expression experiments typically have no more than a few samples for each condition. Current single-cell studies profile hundreds or even thousands of cells in each condition and must therefore be analyzed with algorithms that scale to very large data sets.

A cell-state transition typically involves changes in hundreds or even thousands of genes. How do we identify the “master regulators” of the switch? Inferring regulatory networks from expression data captured with microarrays or RNA-seq is considered extremely hard, but single-cell genomics may make it vastly easier. In theory, variation across cells in the RNA level of a key regulatory gene such as a transcription factor should be reflected in variation in RNA levels of downstream targets. Unfortunately, computational methods for inferring regulatory networks have been hamstrung by two major issues. The first is another manifestation of the curse of dimensionality; because there are so many possible gene–gene interactions, even large-scale experiments lack enough data to reliably predict which genes interact. The second arises due to averaging, which destroys the crucial source of variation that any algorithm needs to accurately reconstruct gene regulatory networks from expression data. Moreover, typical experiments capture hundreds or even thousands of cells under highly controlled, consistent conditions. This constitutes orders of magnitude more data than conventional RNA-seq or microarray studies generate, and thus far more information than any existing network inference algorithm has ever been provided.

Several recent single-cell expression studies have constructed regulatory networks from the data using a variety of machine learning–based methods. For example, Buganim et al. (2012) constructed a Bayesian network to define the hierarchy of transcriptional regulators that drive reprogramming from mouse fibroblasts to the induced pluripotent state. After reconstructing a branched trajectory with diffusion maps, Moignard et al. (2015) synthesized a Boolean regulatory network model to pinpoint the genes driving the trajectory. Both of these studies have some important caveats. Notably, because they probe 48 or 96 different genes, the network might exclude some important regulators. Nevertheless, these studies are likely to be the first of many that aim to realize the long-sought goal of inferring regulatory relationships directly from expression data.

Conclusions and outlook

Single-cell genomics experiments will revolutionize our understanding of gene regulation during development and disease because they access a level of information that has simply never been available. Fundamental questions regarding how cell identity is defined and maintained will be, for the first time, answerable. Robust unsupervised clustering analysis will help catalog cell types and states. Trajectory analysis will track cells as they travel across the landscape of possible states. In time, we may even be able to construct complete regulatory networks and genome-scale dynamical systems of the cell from these measurements. However, much work remains to be done, both in the experimental and computational single-cell domains. Sensitivity of single-cell RNA-seq and other sequencing assays needs further improvement. There are still aspects of genome regulation, such as transcription-factor binding activity, that are not yet measurable in single cells. Many basic statistical analysis questions, such as how to deal with batch effects and technical variability, have not been resolved. Algorithms to identify the critical regulators of cell fate decisions are only just beginning to appear. Despite rapid developments and considerable uncertainty surrounding protocols, instruments, and software, one thing is clear: Single-cell genomics will have a transformative impact on cell biology.