Install packages

This notebook demonstrates the use of command line tools kallisto and bustools. Please use kallisto >= 0.46, whose binary can be downloaded here. The binary of bustools can be found here.

After you download the binary, you should decompress the file (if it is tar.gz) with tar -xzvf file.tar.gz in the bash terminal, and add the directory containing the binary to PATH by export PATH=$PATH:/foo/bar, where /foo/bar is the directory of interest. Then you can directly invoke the binary on the command line as we will do in this notebook.

We will be using the R packages below. BUSpaRse is not yet on CRAN or Bioconductor. For Mac users, see the installation note for BUSpaRse. BUSpaRse will be used to generate the transcript to gene file for bustools and to read output of bustools into R. We will also use Seurat version 3 which is now on CRAN. Recently, Satija lab announced SeuratWrappers, with which we can run RNA velocity directly from Seurat. SeuratWrappers is also GitHub only at present. We need to install velocyto.R, which is GitHub only, to compute and visualize RNA velocity after quantifying spliced and unspliced transcripts.

# Install devtools if it's not already installed
if (!require(devtools)) {
  install.packages("devtools")
}
# Install from GitHub
devtools::install_github("BUStools/BUSpaRse")
devtools::install_github("satijalab/seurat-wrappers")
devtools::install_github("velocyto-team/velocyto.R")

This vignette uses the version of DropletUtils from Bioconductor version 3.9; the version from Bioconductor 3.8 has a different user interface. If you are using a version of R older than 3.6.0 and want to rerun this vignette, then you can adapt the knee plot code to the older version of DropletUtils, or install DropletUtils from GitHub, which I did for this notebook. BSgenome.Mmusculus.UCSC.mm10 and AnnotationHub are also on Bioconductor. Bioconductor packages can be installed as such:

if (!require(BiocManager)) {
  install.packages("BiocManager")
}
BiocManager::install(c("DropletUtils", "BSgenome.Mmusculus.UCSC.mm10", "AnnotationHub"))

The other packages are on CRAN.

library(BUSpaRse)
library(Seurat)
library(SeuratWrappers)
library(BSgenome.Mmusculus.UCSC.mm10)
library(AnnotationHub)
library(zeallot) # For %<-% that unpacks lists in the Python manner
library(DropletUtils)
library(tidyverse)
library(uwot) # For umap
library(GGally) # For ggpairs
library(velocyto.R)
library(SingleR)
library(scales)
library(plotly)
theme_set(theme_bw())

Download data

The dataset we are using is 10x 10k neurons from an E18 mouse (almost 25 GB).

# Download data
if (!file.exists("./data/neuron_10k_fastqs.tar")) {
  download.files("http://s3-us-west-2.amazonaws.com/10x.files/samples/cell-exp/3.0.0/neuron_10k_v3/neuron_10k_v3_fastqs.tar", "./data/neuron_10k_fastqs.tar", method = "wget", quiet = TRUE)
}

Then untar the downloaded file.

cd ./data
tar -xvf ./data/neuron_10k_fastqs.tar

Remove empty droplets

Now we have the spliced and unspliced matrices to be read into R:

c(spliced, unspliced) %<-% read_velocity_output(spliced_dir = "./output/neuron10k_collapse/spliced",
                                                spliced_name = "s",
                                                unspliced_dir = "./output/neuron10k_collapse/unspliced",
                                                unspliced_name = "u")

The %<-% from zeallot unpacks a list of 2 into 2 separate objects in the Python and Matlab manner. How many UMIs are from unspliced transcripts?

sum(unspliced@x) / (sum(unspliced@x) + sum(spliced@x))

[1] 0.5311959

There are more unspliced counts than spliced counts, which has been observed in multiple datasets. In contrast, for velocyto, the unspliced count is usually between 10% and 20% of the sum of spliced and unspliced. Perhaps this is because kallisto | bus counts reads that are partially intronic and partially exonic as unspliced while velocyto throws away many reads (see this GitHub issue).

We expect around 10,000 cells. There are over 10 times more barcodes here, since most barcodes are from empty droplets. The number of genes does not seem too outrageous.

dim(spliced)

[1]   54446 1213262

dim(unspliced)

[1]   54446 1382092

Most barcodes only have 0 or 1 UMIs detected.

tot_count <- Matrix::colSums(spliced)
summary(tot_count)

    Min.  1st Qu.   Median     Mean  3rd Qu.     Max. 
    0.00     1.00     1.00    64.73     2.00 46852.00 

A commonly used method to estimate the number of empty droplets is barcode ranking knee and inflection points, as those are often assumed to represent transition between two components of a distribution. While more sophisticated methods exist (e.g. see emptyDrops in DropletUtils), for simplicity, we will use the barcode ranking method here. However, whichever way we go, we don’t have the ground truth. The spliced matrix is used for filtering, though both matrices have similar inflection points.

bc_rank <- barcodeRanks(spliced)
bc_uns <- barcodeRanks(unspliced)

Here the knee plot is transposed, because this is more generalizable to multi-modal data, such that those with not only RNA-seq but also abundance of cell surface markers. In that case, we can plot number of UMIs on the x axis, number of cell surface protein tags on the y axis, and barcode rank based on both UMI and protein tag counts on the z axis; it makes more sense to make barcode rank the dependent variable. See this blog post by Lior Pachter for a more detailed explanation.

tibble(rank = bc_rank$rank, total = bc_rank$total, matrix = "spliced") %>% 
  bind_rows(tibble(rank = bc_uns$rank, total = bc_uns$total, matrix = "unspliced")) %>% 
  distinct() %>% 
  ggplot(aes(total, rank, color = matrix)) +
  geom_line() +
  geom_vline(xintercept = metadata(bc_rank)$knee, color = "blue", linetype = 2) +
  geom_vline(xintercept = metadata(bc_rank)$inflection, color = "green", linetype = 2) +
  geom_vline(xintercept = metadata(bc_uns)$knee, color = "purple", linetype = 3) +
  geom_vline(xintercept = metadata(bc_uns)$inflection, color = "cyan", linetype = 3) +
  annotate("text", y = c(1000, 1000, 500, 500), 
           x = 1.5 * c(metadata(bc_rank)$knee, metadata(bc_rank)$inflection,
                       metadata(bc_uns)$knee, metadata(bc_uns)$inflectio),
           label = c("knee (s)", "inflection (s)", "knee (u)", "inflection (u)"), 
           color = c("blue", "green", "purple", "cyan")) +
  scale_x_log10() +
  scale_y_log10() +
  labs(y = "Rank", x = "Total UMI counts") +
  theme_bw()

Warning: Transformation introduced infinite values in continuous x-axis

Past versions of knee-1.png

Which inflection point should be used to remove what are supposed to be empty droplets? The one of the spliced matrix or the unspliced matrix?

Actually, spliced and unspliced counts are multimodal data, so why not make one of those promised 3D plots where the barcode rank depends on two variables? The rank (z axis) would now be the number cells with at least x spliced UMIs and y unspliced UMIs. How shall this be computed? The transposed knee plot (or rank-UMI plot) can be thought of as (1 - ECDF(total_UMI))*n_cells. In the ECDF of total UMI counts, the dependent variable is the proportion of cells with at most this number of distinct UMIs. So 1 minus that would mean the proportion of cells with at least this number of distinct UMIs. In the knee plot, the rank is the number of cells with at least this number of distinct UMIs. So dividing by the number of cells, we get 1 - ECDF(total_UMI). Would computing the 2D ECDF be more efficient than this naive approach? There is an R package that can compute bivariate ECDFs called Emcdf, but it uses so much memory that even our server can’t handle. I failed to find implementations of bivariate ECDFs in other languages. There is an algorithm based on range trees that can find multivariate ECDF efficiently.

Before obtaining a more efficient implementation, I used my naive approach that translates this concept into code very literally. Though I used Rcpp, it’s really slow. The trick to make it faster is to only evaluate how many cells have at least x spliced and y unspliced counts at a smaller number of grid points of x and y.

//[[Rcpp::depends(RcppProgress)]]
#include <progress.hpp>
#include <progress_bar.hpp>
#include <Rcpp.h>
using namespace Rcpp;

//[[Rcpp::export]]
NumericMatrix bc_ranks2(NumericVector x, NumericVector y, 
                        NumericVector x_grid, NumericVector y_grid) {
  NumericMatrix out(x_grid.size(), y_grid.size());
  Progress p(x_grid.size(), true);
  for (int i = 0; i < x_grid.size(); i++) {
    checkUserInterrupt();
    for (int j = 0; j < y_grid.size(); j++) {
      out(i,j) = sum((x_grid[i] <= x) & (y_grid[j] <= y));
    }
    p.increment();
  }
  return(out);
}

As most barcodes have a small number of distinct UMIs detected, the grid should be denser for fewer counts. Making the grid in log space achieves this.

# Can only plot barcodes with both spliced and unspliced counts
bcs_inter <- intersect(colnames(spliced), colnames(unspliced))
s <- colSums(spliced[,bcs_inter])
u <- colSums(unspliced[,bcs_inter])
# Grid points
sr <- sort(unique(exp(round(log(s)*100)/100)))
ur <- sort(unique(exp(round(log(u)*100)/100)))

# Run naive approach
bc2 <- bc_ranks2(s, u, sr, ur)

What would the “rank” look like?

# can't turn color to lot scale unless log values are plotted
z_use <- log10(bc2)
z_use[is.infinite(z_use)] <- NA
plot_ly(x = sr, y = ur, z = z_use) %>% add_surface() %>% 
  layout(scene = list(xaxis = list(title = "Total spliced UMIs", type = "log"),
                      yaxis = list(title = "Total unspliced UMIs", type = "log"),
                      zaxis = list(title = "Rank (log10)")))

Looks like it worked. This looks pretty symmetric as the rank-UMI plots for the spliced and unspliced matrices are pretty similar. How can this be used to decide what may be empty droplets? This worths some more thoughts. The surface might also need to be be smoothed for automated thresholding, just like in DropletUtils’s inflection method. For now, for simplicity, the inflection point for the spliced matrix will be used provisionally.

bcs_use <- colnames(spliced)[tot_count > metadata(bc_rank)$inflection]
# Remove genes that aren't detected
tot_genes <- Matrix::rowSums(spliced)
genes_use <- rownames(spliced)[tot_genes > 0]
sf <- spliced[genes_use, bcs_use]
uf <- unspliced[genes_use, bcs_use]

dim(sf)

[1] 24385 12394

Cell type annotation

SingleR uses bulk RNA-seq data of isolated known cell types as a reference to annotate cell types in scRNA-seq datasets. The reference uses gene symbols or names rather than Ensembl IDs.

# Get gene names
gns <- tr2g_EnsDb(edb)[,c("gene", "gene_name")] %>% 
  distinct()

data("mouse.rnaseq")
# Convert from gene symbols to Ensembl gene ID
ref_use <- mouse.rnaseq$data
rownames(ref_use) <- gns$gene[match(rownames(ref_use), gns$gene_name)]
ref_use <- ref_use[!is.na(rownames(ref_use)),]

annot <- SingleR("single", sf, ref_data = ref_use, types = mouse.rnaseq$types)

In order not to have cells not of the neural or glial lineages overshadow velocity visualization, only cells of the neural and glial lineages are kept.

ind <- annot$labels %in% c("NPCs", "Neurons", "OPCs", "Oligodendrocytes", 
                           "qNSCs", "aNSCs", "Astrocytes", "Ependymal")
cells_use <- annot$cell.names[ind]
sf <- sf[, cells_use]
uf <- uf[, cells_use]

Meaning of the acronyms:

• NPCs: Neural progenitor cells
• OPCs: Oligodendrocyte progenitor cells
• qNSCs: Quiescent neural stem cells
• aNSCs: Active neural stem cells

QC

Both the spliced and unspliced matrices are normalized and scaled with SCTransform, which is an alternative to NormalizeData, ScaleData, and FindVariableFeatures.

seu <- CreateSeuratObject(sf, assay = "sf") %>% 
  SCTransform(assay = "sf", new.assay.name = "spliced")
seu[["uf"]] <- CreateAssayObject(uf)
seu <- SCTransform(seu, assay = "uf", new.assay.name = "unspliced")
# Add cell type metadata
seu <- AddMetaData(seu, setNames(annot$labels[ind], cells_use), 
                   col.name = "cell_type")

cols_use <- c("nCount_sf", "nFeature_sf", "nCount_uf", "nFeature_uf")
VlnPlot(seu, cols_use, pt.size = 0.1, ncol = 1, group.by = "cell_type")

Past versions of unnamed-chunk-38-1.png

There’s only one cell labeled ependymal by SingleR. How does number of UMI counts relate to number of genes detected? How does number of UMI counts in the spliced matrix relate to the number of gene detected in the unspliced matrix?

# Helper functions for ggpairs
log10_diagonal <- function(data, mapping, ...) {
  ggally_densityDiag(data, mapping, ...) + scale_x_log10()
}
log10_points <- function(data, mapping, ...) {
  ggally_points(data, mapping, ...) + scale_x_log10() + scale_y_log10()
}

ggpairs(seu@meta.data, columns = cols_use,
        upper = list(continuous = "cor"),
        diag = list(continuous = log10_diagonal),
        lower = list(continuous = wrap(log10_points, alpha = 0.1, size=0.3)),
        progress = FALSE)

Past versions of pairs-1.png

Dimension reduction

When visualizing RNA velocity on reduced dimensions, should the cell embeddings be from the spliced matrix or the unspliced matrix or the sum of both? In my opinion, it makes the most sense to plot RNA velocity over cell embeddings from the spliced matrix. The arrows in RNA velocity visualization stand for where the cell is predicted to be going in the near future. Where does the cell go from? The current state. And the current state is represented by the spliced matrix, while the unspliced matrix represents what is soon to come. Thus all the dimension reduction here will be computed from the spliced matrix.

DefaultAssay(seu) <- "spliced"
seu <- RunPCA(seu, verbose = FALSE, npcs = 70)
ElbowPlot(seu, ndims = 70)

Past versions of elbow-1.png

PCAPlot(seu, group.by = "cell_type", pt.size = 0.5, label = TRUE, repel = TRUE) +
  scale_color_brewer(type = "qual", palette = "Set2")

Warning: Using `as.character()` on a quosure is deprecated as of rlang 0.3.0.
Please use `as_label()` or `as_name()` instead.
This warning is displayed once per session.

Past versions of pca-1.png

seu <- RunTSNE(seu, dims = 1:50, verbose = FALSE)
TSNEPlot(seu, group.by = "cell_type", pt.size = 0.5, label = TRUE, repel = TRUE) +
  scale_color_brewer(type = "qual", palette = "Set2")

Past versions of tsne-1.png

This looks quite similar to the tSNE from gene count matrix of this same dataset, except rotated; see the slingshot notebook

Seurat’s RunUMAP doesn’t work on servers, where this notebook is built, because of some problem with reticulate, which RunUMAP uses to call the Python implementation of UMAP. See this GitHub issue. Thus we will use the R package uwot for UMAP instead.

set.seed(123)
um <- umap(Embeddings(seu, reduction = "pca")[, 1:50], init = "random")
rownames(um) <- colnames(sf)
colnames(um) <- c("UMAP_1", "UMAP_2")
seu[["umap"]] <- CreateDimReducObject(embeddings = um, assay = "spliced", key = "umap_")

UMAPPlot(seu, pt.size = 0.5, group.by = "cell_type", label = TRUE, repel = TRUE) +
  scale_color_brewer(type = "qual", palette = "Set2")

Past versions of umap-1.png

As expected, qNSCs are on the one end , and neurons are on the other. Clustering should partition the big blob of NPCs that SingleR could not further partition due to limitations in the SingleR reference for mouse brains.

seu <- FindNeighbors(seu, verbose = FALSE) %>% 
  FindClusters(resolution = 1, verbose = FALSE) # Louvain
UMAPPlot(seu, pt.size = 0.5)

Past versions of umap_clust-1.png