Slingshot: Lineage Inference for Single-cell Data

knitr::opts_chunk$set(fig.align="center", cache=TRUE,error=FALSE, #stop on error
fig.width=5, fig.height=5, autodep=TRUE, out.width="600px", out.height="600px",
results="markup", echo=TRUE, eval=TRUE)
#knitr::opts_knit$set(stop_on_error = 2L) #really make it stop
#knitr::dep_auto()
options(getClass.msg=FALSE)
graphics:::par(pch = 16, las = 1)
set.seed(12345) ## for reproducibility

library(slingshot, quietly = TRUE)

Introduction

This vignette will demonstrate a full single-cell lineage analysis workflow, with particular emphasis on the processes of lineage reconstruction and pseudotime inference. We will make use of the slingshot package proposed in [@slingshot] and show how it may be applied in a broad range of settings.

The goal of slingshot is to use clusters of cells to uncover global structure and convert this structure into smooth lineages represented by one-dimensional variables, called "pseudotime." We provide tools for learning cluster relationships in an unsupervised or semi-supervised manner and constructing smooth curves representing each lineage, along with visualization methods for each step.

Overview

The minimal input to slingshot is a matrix representing the cells in a reduced-dimensional space and a vector of cluster labels. With these two inputs, we then:

Datasets

We will work with two simulated datasets in this vignette. The first was generated by the splatter package [@splatter], using parameters inferred from the single-cell dataset contained in the HSMMSingleCell package. This dataset is contained in a SingleCellExperiment object [@SingleCellExperiment] and will be used to demonstrate a full "start-to-finish" workflow. It can be loaded via the splatterPath dataset, provided with slingshot. Below, we provide code showing how the dataset was generated.

library(HSMMSingleCell, quietly = TRUE)
library(splatter, quietly = TRUE)

# load data
data("HSMM_expr_matrix")
counts <- round(HSMM_expr_matrix)

# estimate simulation parameters
params <- splatEstimate(counts)

# set path-related parameters
params <- setParams(params, update = list(batchCells = 200, 
                                          path.nonlinearProb = .5, 
                                          de.prob = .6, de.facLoc = 1, 
                                          de.facScale = 2))
library(splatter, quietly = TRUE)
data("splatterParameters")
# simulate
sim <- splatSimulate(params, method = "paths")

The second dataset consists of a $140$ cells $\times$ $2$ dimensions matrix of coordinates, along with cluster labels generated by $k$-means clustering. This dataset represents a bifurcating lineage trajectory and it will allow us to demonstrate some of the additional functionality offered by slingshot.

library(slingshot, quietly = FALSE)
data("slingshotExample")

dim(rd) # data representing cells in a reduced dimensional space
length(cl) # vector of cluster labels

Upstream Analysis

Gene Filtering

To begin our analysis of the single lineage dataset, we need to reduce the dimensionality of our data and filtering out uninformative genes is a typical first step. This will greatly improve the speed of downstream analyses, while keeping the loss of information to a minimum.

For the gene filtering step, we retained any genes robustly expressed in at least enough cells to constitute a cluster, making them potentially interesting cell-type marker genes. We set this minimum cluster size to 10 cells and define a gene as being "robustly expressed" if it has a simulated count of at least 10 reads.

# filter genes down to potential cell-type markers
# at least M (15) reads in at least N (15) cells
geneFilter <- apply(assays(sim)$counts,1,function(x){
    sum(x >= 15) >= 15
})
sim <- sim[geneFilter, ]

Normalization

Another important early step in most RNA-Seq analysis pipelines is the choice of normalization method. This allows us to remove unwanted technical or biological artifacts from the data, such as batch, sequencing depth, cell cycle effects, etc. In practice, it is valuable to compare a variety of normalization techniques and compare them along different evaluation criteria, for which we recommend the scone package [@scone].

Since we are working with simulated data, we have the advantage of knowing what effects are present. In this case, we know that each simulated cell had a unique expected library size and we will apply quantile normalization in order to remove this effect.

FQnorm <- function(counts){
  rk <- apply(counts,2,rank,ties.method='min')
  counts.sort <- apply(counts,2,sort)
  refdist <- apply(counts.sort,1,median)
  norm <- apply(rk,2,function(r){ refdist[r] })
  rownames(norm) <- rownames(counts)
  return(norm)
}
assays(sim)$norm <- FQnorm(assays(sim)$counts)

Dimensionality Reduction

The fundamental assumption of slingshot is that cells which are transcriptionally similar will be close to each other in some reduced-dimensional space. Since we use Euclidean distances in constructing lineages and measuring pseudotime, it is important to have a low-dimensional representation of the data.

There are many methods available for this task and we will intentionally avoid the issue of determining which is the "best" method, as this likely depends on the type of data, method of collection, upstream computational choices, and many other factors. We will demonstrate two methods of dimensionality reduction: principal components analysis (PCA) and diffusion maps (via the destiny package, [@destiny]).

When performing PCA, we do not scale the genes by their variance because we do not believe that all genes are equally informative. We want to find signal in the robustly expressed, highly variable genes, not dampen this signal by forcing equal variance across genes. When plotting, we make sure to set the aspect ratio, so as not to distort the perceived distances.

pca <- prcomp(t(log1p(assays(sim)$norm)), scale. = FALSE)
rd1 <- pca$x[,1:2]

plot(rd1, col = topo.colors(100)[sim$Step], pch=16, asp = 1)
library(destiny, quietly = TRUE)
dm <- DiffusionMap(t(log1p(assays(sim)$norm)))
rd2 <- cbind(DC1 = dm$DC1, DC2 = dm$DC2)

plot(rd2, col = topo.colors(100)[sim$Step], pch=16, asp = 1)

We will add both dimensionality reductions to the SingleCellExperiment object, but continue our analysis focusing on the PCA results.

reducedDims(sim) <- SimpleList(PCA = rd1, DiffMap = rd2)

Clustering Cells

The final input to slingshot is a vector of cluster labels for the cells. If this is not provided, slingshot will treat the data as a single cluster and fit a standard principal curve. However, we recommend clustering the cells even in datasets where only a single lineage is expected, as it allows for the potential discovery of novel branching events.

The clusters identified in this step will be used to determine the global structure of the underlying lineages (that is, their number, when they branch off from one another, and the approximate locations of those branching events). This is different than the typical goal of clustering single-cell data, which is to identify all biologically relevant cell types present in the dataset. For example, when determining global lineage structure, there is no need to distinguish between immature and mature neurons since both cell types will, presumably, fall along the same segment of a lineage.

For our analysis, we implement two clustering methods which similarly assume that Euclidean distance in a low-dimensional space reflect biological differences between cells: Gaussian mixture modeling and $k$-means. The former is implemented in the mclust package [@mclust] and features an automated method for determining the number of clusters based on the Bayesian information criterion (BIC).

library(mclust, quietly = TRUE)
cl1 <- Mclust(rd1)$classification
colData(sim)$GMM <- cl1

library(RColorBrewer)
plot(rd1, col = brewer.pal(9,"Set1")[cl1], pch=16, asp = 1)

While $k$-means does not have a similar functionality, we have shown in [@slingshot] that simultaneous principal curves are quite robust to the choice of $k$, so we select a $k$ of 4 somewhat arbitrarily. If this is too low, we may miss a true branching event and if it is too high or there is an abundance of small clusters, we may begin to see spurious branching events.

cl2 <- kmeans(rd1, centers = 4)$cluster
colData(sim)$kmeans <- cl2

plot(rd1, col = brewer.pal(9,"Set1")[cl2], pch=16, asp = 1)

Using Slingshot

At this point, we have everything we need to run slingshot on our simulated dataset. This is a two-step process composed of \textbf{1.)} identifying the global lineage structure with a cluster-based minimum spanning tree (MST) and \textbf{2.)} fitting simultaneous principal curves to describe each lineage.

These two steps can be run separately with the getLineages and getCurves functions, or together with the wrapper function, slingshot (recommended). We will use the combined method for the analysis of the splatter dataset, but demonstrate the usage of the individual functions later, on the slingshot dataset.

The slingshot wrapper function performs both steps of lineage inference in a single call. The necessary inputs are a reduced dimensional matrix of coordinates and a set of cluster labels. These can be separate objects or, in the case of the splatter data, elements contained in a SingleCellExperiment object.

To run slingshot with the dimensionality reduction produced by PCA and cluster labels identified by Gaussian mixutre modeling, we would do the following:

sce <- slingshot(sim, clusterLabels = 'GMM', reducedDim = 'PCA')

As noted above, if no clustering results are provided, it is assumed that all cells are part of the same cluster and a single curve will be constructed. If no dimensionality reduction is provided, slingshot will use the first element of the list returned by reducedDims.

The output is a SingleCellExperiment object with slingshot results incorporated. Most of the output is added to the metadata in the form of a list and is accessible via metadata(sce)$slingshot. Additionally, all inferred pseudotime variables (one per lineage) are added to the colData. To extract all slingshot results in a single object, we can use the SlingshotDataSet function. This can be useful for visualization, as several plotting methods for SlingshotDataSet objects are included with the package. Below, we visuzalize the inferred lineage for the splatter data with points colored by pseudotime.

summary(sce$slingPseudotime_1)

colors <- colorRampPalette(brewer.pal(11,'Spectral')[-6])(100)
plot(reducedDims(sce)$PCA, col = colors[cut(sce$slingPseudotime_1,breaks=100)], pch=16, asp = 1)
lines(SlingshotDataSet(sce), lwd=2)

We can also see how the lineage structure was intially estimated by the cluster-based minimum spanning tree by using the type argument.

plot(reducedDims(sce)$PCA, col = brewer.pal(9,'Set1')[sce$GMM], pch=16, asp = 1)
lines(SlingshotDataSet(sce), lwd=2, type = 'lineages')

Downstream Analysis

Identifying temporally expressed genes

After running slingshot, an interesting next step may be to find genes that change their expression over the course of development. We demonstrate one possible method for this type of analysis on the $1,000$ most variable genes.

We will regress each gene on the pseudotime variable we have generated, using a general additive model (GAM). This allows us to detect non-linear patterns in gene expression.

require(gam)
t <- sce$slingPseudotime_1

# for time, only look at the 1,000 most variable genes
Y <- log1p(assays(sim)$norm)
var1K <- names(sort(apply(Y,1,var),decreasing = TRUE))[1:1000]
Y <- Y[var1K,]

# fit a GAM with a loess term for pseudotime
gam.pval <- apply(Y,1,function(z){
  d <- data.frame(z=z, t=t)
  tmp <- gam(z ~ lo(t), data=d)
  p <- summary(tmp)[4][[1]][1,5]
  p
})

We can then pick out the top genes based on p-value and visualize their expression over developmental time with a heatmap.

require(clusterExperiment)
topgenes <- names(sort(gam.pval, decreasing = FALSE))[1:100]
heatdata <- assays(sim)$norm[rownames(assays(sim)$norm) %in% topgenes, 
                        order(t, na.last = NA)]
heatclus <- sce$GMM[order(t, na.last = NA)]
ce <- ClusterExperiment(heatdata, heatclus, transformation = log1p)
plotHeatmap(ce, clusterSamplesData = "orderSamplesValue",visualizeData = 'transformed')

Detailed Slingshot Functionality

Here, we provide further details and highlight some additional functionality of the slingshot package. We will use the included slingshotExample dataset for illustrative purposes. This dataset was designed to represent cells in a low dimensional space and comes with a set of cluster labels generated by $k$-means clustering. Rather than constructing a full SingleCellExperiment object, which requires gene-level data, we will use the low-dimensional matrix of coordinates directly and provide the cluster labels as an additional argument.

Identifying global lineage structure

The getLineages function takes as input an n $\times$ p matrix and a vector of clustering results of length n. It maps connections between adjacent clusters using a minimum spanning tree (MST) and identifies paths through these connections that represent lineages. The output of this function is a SlingshotDataSet containing the inputs as well as the inferred MST (represented by an adjacency matrix) and lineages (ordered vectors of cluster names).

This analysis can be performed in an entirely unsupervised manner or in a semi-supervised manner by specifying known initial and terminal point clusters. If we do not specify a starting point, slingshot selects one based on parsimony, maximizing the number of clusters shared between lineages before a split. If there are no splits or multiple clusters produce the same parsimony score, the starting cluster is chosen arbitrarily.

In the case of our simulated data, slingshot selects Cluster 1 as the starting cluster. However, we generally recommend the specification of an initial cluster based on prior knowledge (either time of sample collection or established gene markers). This specification will have no effect on how the MST is constructed, but it will impact how branching curves are constructed.

lin1 <- getLineages(rd, cl, start.clus = '1')
lin1
plot(rd, col = brewer.pal(9,"Set1")[cl], asp = 1, pch = 16)
lines(lin1, lwd = 3)

At this step, slingshot also allows for the specification of known endpoints. Clusters which are specified as terminal cell states will be constrained to have only one connection when the MST is constructed (ie., they must be leaf nodes). This constraint could potentially impact how other parts of the tree are drawn, as shown in the next example where we specify Cluster 3 as an endpoint.

lin2 <- getLineages(rd, cl, start.clus= '1', end.clus = '3')

plot(rd, col = brewer.pal(9,"Set1")[cl], asp = 1, pch = 16)
lines(lin2, lwd = 3, show.constraints = TRUE)

This type of supervision can be useful for ensuring that results are consistent with prior biological knowledge. Specifically, it can prevent known terminal cell fates from being categorized as transitory states.

There are a few additional arguments we could have passed to getLineages for greater control:

After constructing the MST, getLineages identifies paths through the tree to designate as lineages. At this stage, a lineage will consist of an ordered set of cluster names, starting with the root cluster and ending with a leaf. The output of getLineages is a SlingshotDataSet which holds all of the inputs as well as a list of lineages and some additional information on how they were constructed. This object is then used as the input to the getCurves function.

Constructing smooth curves and ordering cells

In order to model development along these various lineages, we will construct smooth curves with the function getCurves. Using smooth curves based on all the cells eliminates the problem of cells projecting onto vertices of piece-wise linear trajectories and makes slingshot more robust to noise in the clustering results.

In order to construct smooth lineages, getCurves follows an iterative process similar to that of principal curves presented in [@princurve]. When there is only a single lineage, the resulting curve is simply the principal curve through the center of the data, with one adjustment: the initial curve is constructed with the linear connections between cluster centers rather than the first prinicpal component of the data. This adjustment adds stability and typically hastens the algorithm's convergence.

When there are two or more lineages, we add an additional step to the algorithm: averaging curves near shared cells. Both lineages should agree fairly well on cells that have yet to differentiate, so at each iteration we average the curves in the neighborhood of these cells. This increases the stability of the algorithm and produces smooth branching lineages.

crv1 <- getCurves(lin1)
crv1
plot(rd, col = brewer.pal(9,"Set1")[cl], asp = 1, pch = 16)
lines(crv1, lwd = 3)

The output of getCurves is an updated slingshotDataSet which now contains the simultaneous principal curves and additional information on how they were fit. The pseudotime function extracts a cells-by-lineages matrix of pseudotime values for each cell along each lineage, with NA values for cells that were not assigned to a particular lineage. The curveWeights function extracts a similar matrix of weights assigning each cell to one or more lineages.

The curves objects can be accessed using the curves function, which will return a list of principal_curve objects. These objects consist of the following slots:

Session Info

sessionInfo()

References



Try the slingshot package in your browser

Any scripts or data that you put into this service are public.

slingshot documentation built on Nov. 1, 2018, 3 a.m.