Diffusion maps are spectral method for non-linear dimension reduction introduced by Coifman et al. (2005). Diffusion maps are based on a distance metric (diffusion distance) which is conceptually relevant to how differentiating cells follow noisy diffusion-like dynamics, moving from a pluripotent state towards more differentiated states.

The R package destiny implements the formulation of diffusion maps presented in Haghverdi, Buettner, and Theis (2015) which is especially suited for analyzing single-cell gene expression data from time-course experiments. It implicitly arranges cells along their developmental path, with bifurcations where differentiation events occur.

In particular we follow Haghverdi, Buettner, and Theis (2015) and present an implementation of diffusion maps in R that is less affected by sampling density heterogeneities and thus capable of identifying both abundant and rare cell populations. In addition, destiny implements complex noise models reflecting zero-inflation/censoring due to drop-out events in single-cell qPCR data and allows for missing values. Finally, we further improve on the implementation from Haghverdi, Buettner, and Theis (2015), and implement a nearest neighbour approximation capable of handling very large data sets of up to 300.000 cells.

For those familiar with R, and data preprocessing, we recommend the section Plotting.

Preprocessing of single qPCR data

As an example, we present in the following the preprocessing of data from Guo et al. (2010). This dataset was produced by the Biomark RT-qPCR system and contains Ct values for 48 genes of 442 mouse embryonic stem cells at 7 different developmental time points, from the zygote to blastocyst.

Starting at the totipotent 1-cell stage, cells transition smoothly in the transcriptional landscape towards either the trophoectoderm lineage or the inner cell mass. Subsequently, cells transition from the inner cell mass either towards the endoderm or epiblast lineage. This smooth transition from one developmental state to another, including two bifurcation events, is reflected in the expression profiles of the cells and can be visualized using destiny.

Import

Downloading the table S4 from the publication website will give you a spreadsheet “mmc4.xls”, from which the data can be loaded:

library(readxl)
raw_ct <- read_xls('mmc4.xls', 'Sheet1')

raw_ct[1:9, 1:9]  #preview of a few rows and columns

## # A tibble: 9 × 9
##   Cell   Actb  Ahcy  Aqp3 Atp12a  Bmp4  Cdx2 Creb312 Cebpa
##   <chr> <dbl> <dbl> <dbl>  <dbl> <dbl> <dbl>   <dbl> <dbl>
## 1 1C 1   14.0  19.3  23.9     28    28  21.3    20.8  28  
## 2 1C 2   13.7  18.6  28       28    28  23.4    20.9  28  
## 3 1C 3   13.4  18.2  26.2     28    28  22.9    19.6  28  
## 4 1C 4   13.7  18.6  28       28    28  23.3    20.7  28  
## 5 1C 5   13.5  18.6  24.2     28    28  24.2    21.8  23.7
## 6 1C 6   12.9  17.4  25.5     28    28  21.9    21.3  28  
## 7 1C 7   13.0  17.4  23.9     28    28  22.7    21.1  28  
## 8 1C 8   12.8  18.4  23.7     28    28  24.1    19.8  28  
## 9 1C 9   13.3  18.3  28       28    28  21.9    21.2  28

The value 28 is the assumed background expression of 28 cycle times.

In order to easily clean and normalize the data without mangling the annotations, we convert the data.frame into a Bioconductor ExpressionSet using the as.ExpressionSet function from the destiny package, and assign it to the name ct:

library(destiny)
library(Biobase)

# as.data.frame because ExpressionSets
# do not support readxl’s “tibbles”
ct <- as.ExpressionSet(as.data.frame(raw_ct))
ct

## ExpressionSet (storageMode: lockedEnvironment)
## assayData: 48 features, 442 samples 
##   element names: exprs 
## protocolData: none
## phenoData
##   sampleNames: 1 2 ... 442 (442 total)
##   varLabels: Cell
##   varMetadata: labelDescription
## featureData: none
## experimentData: use 'experimentData(object)'
## Annotation:

The advantage of ExpressionSet over data.frame is that tasks like normalizing the expressions are both faster and do not accidentally destroy annotations by applying “normalization” to columns that are not expressions. The approach of handling a separate expression matrix and annotation data.frame requires you to be careful when adding or removing samples on both variables, while ExpressionSet does it internally for you.

The object internally stores an expression matrix of features × samples, retrievable using exprs(ct), and an annotation data.frame of samples × annotations as phenoData(ct). Annotations can be accessed directly via ct$column and ct[['column']]. Note that the expression matrix is transposed compared to the usual samples × features data.frame.

Data cleaning

We remove all cells that have a value bigger than the background expression, indicating data points not available (NA). Also we remove cells from the 1 cell stage embryos, since they were treated systematically different (Guo et al. 2010).

For this, we add an annotation column containing the embryonic cell stage for each sample by extracting the number of cells from the “Cell” annotation column:

num_cells <- gsub('^(\\d+)C.*$', '\\1', ct$Cell)
ct$num_cells <- as.integer(num_cells)

We then use the new annotation column to create two filters:

# cells from 2+ cell embryos
have_duplications <- ct$num_cells > 1
# cells with values <= 28
normal_vals <- apply(exprs(ct), 2, function(smp) all(smp <= 28))

We can use the combination of both filters to exclude both non-divided cells and such containing Ct values higher than the baseline, and store the result as cleaned_ct:

cleaned_ct <- ct[, have_duplications & normal_vals]

Normalization

Finally we follow Guo et al. (2010) and normalize each cell using the endogenous controls Actb and Gapdh by subtracting their average expression for each cell. Note that it is not clear how to normalise sc-qPCR data as the expression of housekeeping genes is also stochastic. Consequently, if such housekeeping normalisation is performed, it is crucial to take the mean of several genes.

housekeepers <- c('Actb', 'Gapdh')  # houskeeper gene names

normalizations <- colMeans(exprs(cleaned_ct)[housekeepers, ])

guo_norm <- cleaned_ct
exprs(guo_norm) <- exprs(guo_norm) - normalizations

The resulting ExpressionSet contains the normalized Ct values of all cells retained after cleaning.

Plotting

The data necessary to create a diffusion map with our package is a a cell×gene matrix or data.frame, or alternatively an ExpressionSet (which has a gene×cell exprs matrix). In order to create a DiffusionMap object, you just need to supply one of those formats as first parameter to the DiffusionMap function. In the case of a data.frame, each floating point column is interpreted as expression levels, and columns of different type (e.g. factor, character or integer) are assumed to be annotations and ignored. Note that single-cell RNA-seq count data should be transformed using a variance-stabilizing transformation (e.g. log or rlog); the Ct scale for qPCR data is logarithmic already (an increase in 1 Ct corresponds to a doubling of transcripts).

In order to create a diffusion map to plot, you have to call DiffusionMap, optionally with parameters. If the number of cells is small enough (< ~1000), you do not need to specify approximations like k (for k nearest neighbors).

If you started reading here, execute data(guo.norm) to load the dataset that was created in the previous section.

library(destiny)
#If you start here, run: data(guo_norm)
dm <- DiffusionMap(guo_norm)

Simply calling plot on the resulting object dif will visualize the diffusion components:

plot(dm)

The diffusion map nicely illustrates a branching during the first days of embryonic development.

The annotation column containing the cell stage can be used to annotate our diffusion map. Using the annotation as col parameter will automatically color the map using the current R color palette. Use palette(colors) to configure it globally or plot(..., pal = colors) for one plot.

palette(cube_helix(6)) # configure color palette

plot(dm,
    pch = 20,             # pch for prettier points
    col_by = 'num_cells', # or “col” with a vector or one color
    legend_main = 'Cell stage')

Three branches appear in the map, with a bifurcation occurring the 16 cell stage and the 32 cell stage. The diffusion map is able to arrange cells according to their expression profile: Cells that divided the most and the least appear at the tips of the different branches.

In order to display a 2D plot we supply a vector containing two diffusion component numbers (here 1 & 2) as second argument.

plot(dm, 1:2,
    col_by = 'num_cells',
    legend_main = 'Cell stage')

Alternative visualizations

Diffusion maps consist of eigenvectors called Diffusion Components (DCs) and corresponding eigenvalues. Per default, the first 20 are returned.

You are also able to use packages for interative plots like rgl in a similar fashion, by directly subsetting the DCs using eigenvectors(dif):

if (require(rgl)) {
    plot3d(
        eigenvectors(dm)[, 1:3],
        col = log2(guo_norm$num_cells),
        type = 's', radius = .01)
    view3d(theta = 10, phi = 30, zoom = .8)
    # now use your mouse to rotate the plot in the window
    close3d()
}

## Loading required package: rgl

For the popular ggplot2 package, there is built in support in the form of a fortify.DiffusionMap method, which allows to use DiffusionMap objects as data parameter in the ggplot and qplot functions. But the regular plot function already returns a ggplot object anyway.

library(ggplot2)
ggplot(dm, aes(DC1, DC2, colour = factor(num_cells))) +
    geom_point() + scale_color_cube_helix()

As aesthetics, all diffusion components, gene expressions, and annotations are available. If you plan to make many plots, create a data.frame first by using as.data.frame(dif) or fortify(dif), assign it to a variable name, and use it for plotting.

destiny main vignette: Start here!

Preprocessing of single qPCR data

Import

Data cleaning

Normalization

Plotting

Alternative visualizations

Parameter selection

Dimensions `dims`

destiny main vignette: Start here!

Preprocessing of single qPCR data

Import

Data cleaning

Normalization

Plotting

Alternative visualizations

Parameter selection

Dimensions dims

Dimensions `dims`