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.

# Install devtools if it's not already installed
if (!require(devtools)) {
  install.packages("devtools")
}
# Install from GitHub
devtools::install_github("BUStools/BUSpaRse")

The package SingleR, which will be used for cell type inference, is only available on GitHub:

devtools::install_github("dviraran/SingleR")

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. The package monocle should also be installed from Bioconductor:

if (!require(BiocManager)) {
  install.packages("BiocManager")
}
BiocManager::install(c("DropletUtils", "monocle"))

The other R packages below are on CRAN.

library(BUSpaRse)
library(DropletUtils)
library(monocle)
library(SingleR)
library(Matrix)
library(tidyverse)
library(ggsci)
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)
}

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

Build the kallisto index

Here we use kallisto to pseudoalign the reads to the transcriptome and then to create the bus file to be converted to a sparse matrix. The first step is to build an index of the mouse transcriptome. The transcriptome downloaded here are from Ensembl version 94, released in October 2018.

# Mouse transcriptome
if (!file.exists("./data/mm_cdna.fa.gz")) {
  download.file("ftp://ftp.ensembl.org/pub/release-94/fasta/mus_musculus/cdna/Mus_musculus.GRCm38.cdna.all.fa.gz", "./data/mm_cdna.fa.gz", method = "wget", quiet = TRUE)
}

kallisto index -i ./output/mm_tr_index.idx ./data/mm_cdna.fa.gz

cd ./data/neuron_10k_v3_fastqs
kallisto bus -i ../../output/mm_tr_index.idx -o ../../output/neuron10k -x 10xv3 -t8 \
neuron_10k_v3_S1_L002_R1_001.fastq.gz neuron_10k_v3_S1_L002_R2_001.fastq.gz \
neuron_10k_v3_S1_L001_R1_001.fastq.gz neuron_10k_v3_S1_L001_R2_001.fastq.gz

Run bustools

tr2g <- transcript2gene(fasta_file = "./data/mm_cdna.fa.gz", 
                        kallisto_out_path = "./output/neuron10k",
                        verbose = FALSE)

head(tr2g)

#>              transcript                 gene gene_name
#> 1: ENSMUST00000177564.1 ENSMUSG00000096176.1     Trdd2
#> 2: ENSMUST00000196221.1 ENSMUSG00000096749.2     Trdd1
#> 3: ENSMUST00000179664.1 ENSMUSG00000096749.2     Trdd1
#> 4: ENSMUST00000178862.1 ENSMUSG00000094569.1     Trbd2
#> 5: ENSMUST00000178537.1 ENSMUSG00000095668.1     Trbd1
#> 6: ENSMUST00000179520.1 ENSMUSG00000094028.1   Ighd4-1

# Write tr2g to format required by bustools
save_tr2g_bustools(tr2g, file_save = "./output/tr2g_mm.tsv")

cp ~/cellranger-3.0.2/cellranger-cs/3.0.2/lib/python/cellranger/barcodes/3M-february-2018.txt.gz \
./data/whitelist_v3.txt.gz
gunzip ./data/whitelist_v3.txt.gz

mkdir ./tmp
bustools correct -w ./data/whitelist_v3.txt -p ./output/neuron10k/output.bus | \
bustools sort -T tmp/ -t4 -p - | \
bustools count -o ./output/neuron10k/genes -g ./output/tr2g_mm.tsv \
-e ./output/neuron10k/matrix.ec -t ./output/neuron10k/transcripts.txt --genecounts -
rm -r ./tmp

Remove empty droplets

dim(res_mat)

#> [1]   36121 1410905

The number of genes seems reasonable. The number of barcodes is way larger than the expected ~10k.

tot_counts <- Matrix::colSums(res_mat)
summary(tot_counts)

#>     Min.  1st Qu.   Median     Mean  3rd Qu.     Max. 
#>     0.00     1.00     1.00    64.84     2.00 51473.00

The vast majority of “cells” have only no or just a few UMI detected. Those are empty droplets. 10x claims to have cell capture rate of up to 65%, but in practice, depending on how many cells are in fact loaded, the rate can be much lower. 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.

# Compute barcode rank
bc_rank <- barcodeRanks(res_mat, lower = 1000)

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(total = bc_rank$total, rank = bc_rank$rank) %>% 
  distinct() %>% 
  ggplot(aes(total, rank)) +
  geom_line() +
  geom_vline(xintercept = metadata(bc_rank)$knee, color = "blue", linetype = 2) +
  geom_vline(xintercept = metadata(bc_rank)$inflection, color = "green", linetype = 2) +
  annotate("text", y = 1000, x = 1.5 * c(metadata(bc_rank)$knee, metadata(bc_rank)$inflection),
           label = c("knee", "inflection"), color = c("blue", "green")) +
  scale_x_log10() +
  scale_y_log10() +
  labs(y = "Barcode rank", x = "Total UMI count")

#> Warning: Transformation introduced infinite values in continuous x-axis

# Remove genes that are not detected and empty droplets
res_mat <- res_mat[Matrix::rowSums(res_mat) > 0, tot_counts > metadata(bc_rank)$inflection]
dim(res_mat)

#> [1] 23342 11031

Now the number of cells is closer to expectation.

Cell type inference

Monocle 2 only infers one trajectory for the entire dataset, so non-neuronal cells like endothelial cells and erythrocytes may be mistaken as highly differentiated cells from the neuronal lineage. So we will remove non-neuronal cell types. Cell types are also helpful to orient the trajectory; neuronal progenitor cells must come before neurons. Here cell type inference is done with SingleR, which compares gene expression profiles of individual cells to bulk RNA-seq data of purified known cell types.

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

Then SingleR will assign each cell a label based on Spearman correlation with known cell types from bulk RNA-seq. These are the labels we want:

• NPCs: Neuronal progenitor cells
• aNSCs: Active neuronal stem cells
• qNSCs: Quiescent neuronal stem cells
• Neurons

annots <- SingleR("single", res_mat, ref_data = ref_use, types = mouse.rnaseq$types)
inds <- annots$labels %in% c("NPCs", "aNSCs", "qNSCs", "Neurons")
# Only keep these 4 cell types
cells_use <- annots$cell.names[inds]
mat_neu <- res_mat[,cells_use]

QC

df <- data.frame(nGene = Matrix::colSums(mat_neu > 0),
                 nUMIs = Matrix::colSums(mat_neu))

df %>% 
  gather(key = "variable") %>% 
  ggplot(aes(variable, value)) +
  geom_violin() +
  geom_jitter(size = 0.1, alpha = 0.1, width = 0.3) +
  scale_y_log10()

ggplot(df, aes(nUMIs, nGene)) +
  geom_point(size = 0.5, alpha = 0.3) +
  scale_x_log10() +
  scale_y_log10()