---
title: "Slingshot: Lineage Inference for Single-Cell Data"
author: "Kelly Street"
date: "`r Sys.Date()`"
bibliography: bibFile.bib
output: 
  BiocStyle::html_document:
    toc: true
vignette: >
  %\VignetteEncoding{UTF-8}
---

<!--
%\VignetteEngine{knitr::rmarkdown}
%\VignetteIndexEntry{Slingshot}
-->

```{r options, results="hide", include=FALSE, cache=FALSE, message=FALSE}
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(SingleCellExperiment)
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:

* Identify the global lineage structure by constructing an minimum spanning tree
(MST) on the clusters, with the `getLineages` function.
* Construct smooth lineages and infer pseudotime variables by fitting 
simultaneous principal curves with the `getCurves` function.
* Assess the output of each step with built-in visualization tools.


## Datasets

We will work with two simulated datasets in this vignette. The first (referred
to as the "single-trajectory" dataset) is generated below and designed to
represent a single lineage in which one third of the genes are associated with
the transition. This dataset will be contained in a `SingleCellExperiment`
object [@SingleCellExperiment] and will be used to demonstrate a full
"start-to-finish" workflow.

```{r dataSetup_sim}
# generate synthetic count data representing a single lineage
means <- rbind(
    # non-DE genes
    matrix(rep(rep(c(0.1,0.5,1,2,3), each = 300),100),
        ncol = 300, byrow = TRUE),
    # early deactivation
    matrix(rep(exp(atan( ((300:1)-200)/50 )),50), ncol = 300, byrow = TRUE),
    # late deactivation
    matrix(rep(exp(atan( ((300:1)-100)/50 )),50), ncol = 300, byrow = TRUE),
    # early activation
    matrix(rep(exp(atan( ((1:300)-100)/50 )),50), ncol = 300, byrow = TRUE),
    # late activation
    matrix(rep(exp(atan( ((1:300)-200)/50 )),50), ncol = 300, byrow = TRUE),
    # transient
    matrix(rep(exp(atan( c((1:100)/33, rep(3,100), (100:1)/33) )),50), 
        ncol = 300, byrow = TRUE)
)
counts <- apply(means,2,function(cell_means){
    total <- rnbinom(1, mu = 7500, size = 4)
    rmultinom(1, total, cell_means)
})
rownames(counts) <- paste0('G',1:750)
colnames(counts) <- paste0('c',1:300)
sim <- SingleCellExperiment(assays = List(counts = counts))
```

The second dataset (the "bifurcating" dataset) consists of a matrix of
coordinates (as if obtained by PCA, ICA, diffusion maps, etc.) along with
cluster labels generated by $k$-means clustering. This dataset represents a
bifurcating trajectory and it will allow us to demonstrate some of the
additional functionality offered by `slingshot`.

```{r data_sling}
library(slingshot, quietly = FALSE)
data("slingshotExample")
rd <- slingshotExample$rd
cl <- slingshotExample$cl

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 3
reads.

```{r genefilt}
# 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 >= 3) >= 10
})
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]. We also note that the order of these
steps may change depending on the choice of method. ZINB-WaVE [@zinbwave]
performs dimensionality reduction while accounting for technical variables and
MNN [@mnn] corrects for batch effects after dimensionality reduction.

Since we are working with simulated data, we do not need to worry about batch
effects or other potential confounders. Hence, we will proceed with full
quantile normalization, a well-established method which forces each cell to have
the same distribution of expression values.

```{r norm}
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.

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

plot(rd1, col = rgb(0,0,0,.5), pch=16, asp = 1)
```
```{r dm, cache=TRUE}
library(destiny, quietly = TRUE)
dm <- DiffusionMap(t(log1p(assays(sim)$norm)))
rd2 <- cbind(DC1 = dm$DC1, DC2 = dm$DC2)

plot(rd2, col = rgb(0,0,0,.5), pch=16, asp = 1)
```

We will add both dimensionality reductions to the `SingleCellExperiment` object,
but continue our analysis focusing on the PCA results.
```{r add_RDs, cache=TRUE}
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).

```{r clustering_mclust}
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.

```{r clustering}
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 wrapper function for the analysis of the single-trajectory dataset,
but demonstrate the usage of the individual functions later, on the bifurcating
dataset.

The `slingshot` wrapper function performs both steps of trajectory 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 single-trajectory 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:

```{r sling_sce}
sim <- 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 single-trajectory data with points colored by
pseudotime.

```{r plot_curve_1}
summary(sim$slingPseudotime_1)

colors <- colorRampPalette(brewer.pal(11,'Spectral')[-6])(100)
plotcol <- colors[cut(sim$slingPseudotime_1, breaks=100)]

plot(reducedDims(sim)$PCA, col = plotcol, pch=16, asp = 1)
lines(SlingshotDataSet(sim), lwd=2, col='black')
```

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

```{r plot_curve_2}
plot(reducedDims(sim)$PCA, col = brewer.pal(9,'Set1')[sim$GMM], pch=16, asp = 1)
lines(SlingshotDataSet(sim), lwd=2, type = 'lineages', col = 'black')
```

# 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 $100$ 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.

```{r fitgam, cache = TRUE}
require(gam)
t <- sim$slingPseudotime_1

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

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

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

```{r heatmaps}
topgenes <- names(sort(gam.pval, decreasing = FALSE))[1:100]
heatdata <- assays(sim)$norm[topgenes, order(t, na.last = NA)]
heatclus <- sim$GMM[order(t, na.last = NA)]

heatmap(log1p(heatdata), Colv = NA,
        ColSideColors = brewer.pal(9,"Set1")[heatclus])
```

# 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.

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

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.

```{r lines_sup_end}
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, col = 'black', 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:

* `dist.fun` is a function for computing distances between clusters. The default
is squared distance between cluster centers normalized by their joint covariance
matrix. 
* `omega` is a granularity parameter, allowing the user to set an upper
limit on connection distances. This can be useful for identifying outlier
clusters which are not a part of any lineage.

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.

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

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:

* `s`: the matrix of points that make up the curve. These correspond to the 
orthogonal projections of the data points. 
* `ord`: indices which can be used to
put the cells along a curve in order based their projections. 
* `lambda`: arclengths along the curve from the beginning to each cell's
projection. A matrix of these values is returned by the `slingPseudotime`
function.
* `dist_ind`: the squared distances between data points and their projections
onto the curve.
* `dist`: the sum of the squared projection distances.
* `w`: the vector of weights along this lineage. Cells that were unambiguously
assigned to this lineage will have a weight of `1`, while cells assigned to
other lineages will have a weight of `0`. It is possible for cells to have
weights of `1` (or very close to 1) for multiple lineages, if they are 
positioned before a branching event. A matrix of these values is returned by the
`slingCurveWeights` function.


## Running Slingshot on large datasets

For large datasets, we recommend using the `approx_points` argument with
`slingshot` (or similarly with `getCurves`). While the MST construction operates
on clusters, the process of iteratively projecting all points onto one or more
curves may become computationally impractical as the size of the dataset grows.
By default, `slingshot` constructs curves consisting of $n$ points, where $n$ is
the number of cells, meaning that the computational complexity of the projection
step is proportional to $n^2$. In the presence of branching lineages, these
"dense" curves also affect the complexity of curve averaging and shrinkage.

Setting `approx_points` allows the user to specify the resolution (the number of
unique points) of the curves, which can greatly speed up the computation time
with minimal loss of accuracy. We recommend a value of 100-200. Note that
restricting the number of unique points along the curve does not impose a
similar limit on the number of unique pseudotime values, as demonstrated below
(even with an unrealistically low value of 5 for `approx_points`, we still see a
full color gradient from the pseudotime values):

```{r sling_approxpoints}
sim5 <- slingshot(sim, clusterLabels = 'GMM', reducedDim = 'PCA',
                   approx_points = 5)

colors <- colorRampPalette(brewer.pal(11,'Spectral')[-6])(100)
plotcol <- colors[cut(sim5$slingPseudotime_1, breaks=100)]

plot(reducedDims(sim5)$PCA, col = plotcol, pch=16, asp = 1)
lines(SlingshotDataSet(sim5), lwd=2, col='black')
```


# Session Info

```{r session}
sessionInfo()
```

# References

