Last updated: 2019-02-28

Checks: 5 1

Knit directory: BUSpaRse_notebooks/

This reproducible R Markdown analysis was created with workflowr (version 1.2.0). The Report tab describes the reproducibility checks that were applied when the results were created. The Past versions tab lists the development history.


The R Markdown file has unstaged changes. To know which version of the R Markdown file created these results, you’ll want to first commit it to the Git repo. If you’re still working on the analysis, you can ignore this warning. When you’re finished, you can run wflow_publish to commit the R Markdown file and build the HTML.

Great job! The global environment was empty. Objects defined in the global environment can affect the analysis in your R Markdown file in unknown ways. For reproduciblity it’s best to always run the code in an empty environment.

The command set.seed(20181214) was run prior to running the code in the R Markdown file. Setting a seed ensures that any results that rely on randomness, e.g. subsampling or permutations, are reproducible.

Great job! Recording the operating system, R version, and package versions is critical for reproducibility.

Nice! There were no cached chunks for this analysis, so you can be confident that you successfully produced the results during this run.

Great! You are using Git for version control. Tracking code development and connecting the code version to the results is critical for reproducibility. The version displayed above was the version of the Git repository at the time these results were generated.

Note that you need to be careful to ensure that all relevant files for the analysis have been committed to Git prior to generating the results (you can use wflow_publish or wflow_git_commit). workflowr only checks the R Markdown file, but you know if there are other scripts or data files that it depends on. Below is the status of the Git repository when the results were generated:


Ignored files:
    Ignored:    .Rhistory
    Ignored:    .Rproj.user/
    Ignored:    BUSpaRse_notebooks.Rproj
    Ignored:    data/fastqs/
    Ignored:    data/hgmm_100_fastqs.tar
    Ignored:    data/hgmm_1k_v3_fastqs.tar
    Ignored:    data/hgmm_1k_v3_fastqs/
    Ignored:    data/hs_cdna.fa.gz
    Ignored:    data/mm_cdna.fa.gz
    Ignored:    data/whitelist_v2.txt
    Ignored:    data/whitelist_v3.txt.gz
    Ignored:    output/hs_mm_tr_index.idx
    Ignored:    output/out_hgmm1k/
    Ignored:    output/out_hgmm1k_v3/

Untracked files:
    Untracked:  Rplots.pdf

Unstaged changes:
    Modified:   analysis/10xv2.Rmd
    Modified:   analysis/10xv3.Rmd

Note that any generated files, e.g. HTML, png, CSS, etc., are not included in this status report because it is ok for generated content to have uncommitted changes.


These are the previous versions of the R Markdown and HTML files. If you’ve configured a remote Git repository (see ?wflow_git_remote), click on the hyperlinks in the table below to view them.

File Version Author Date Message
html c465ce1 Lambda Moses 2019-02-14 Build site
Rmd f22beac Lambda Moses 2019-02-14 Added CellRanger whiltelist clarification
Rmd 000215e Lambda Moses 2019-02-14 Clarified git cloning this repo and resolved swapped code chunks for output.sorted.bus
html 71513e5 Lambda Moses 2019-02-14 Good site with figures and reproducibility metrics
Rmd 95f2951 Lambda Moses 2019-02-14 Added head of bus file
html 7184968 Lambda Moses 2019-02-14 Build site.
Rmd 8cd81e0 Lambda Moses 2019-02-14 Finally saved the figures
html 8cd81e0 Lambda Moses 2019-02-14 Finally saved the figures
html 45c5d4a Lambda Moses 2019-02-14 Build site.
Rmd 92aa915 Lambda Moses 2019-02-14 Named chunks with images
html 87d15f5 Lambda Moses 2019-02-14 Build site.
Rmd 8d9fe9a Lambda Moses 2019-02-14 More detailed explanations for kallisto bus workshop
html b7b21a0 Lambda Moses 2019-02-02 Build site.
html a0becdf Lambda Moses 2018-12-19 Build site.
html 93c1053 Lambda Moses 2018-12-19 Build site.
html 85a5770 Lambda Moses 2018-12-14 Build site.
Rmd ca1d6ce Lambda Moses 2018-12-14 Clean up
html 6276894 Lambda Moses 2018-12-14 Build site.
html db70187 Lambda Moses 2018-12-14 Build site.
Rmd b67fc92 Lambda Moses 2018-12-14 Added 10xv3 notebook and further elaboration to 10xv2 notebook
html fff442d Lambda Moses 2018-12-14 Build site.
Rmd 695e202 Lambda Moses 2018-12-14 Changed name to BUSpaRse
html 09a56c1 Lambda Moses 2018-12-14 Build site.
Rmd d288e19 Lambda Moses 2018-12-14 Cache rather than skip evaluation
html 30c9aa3 Lambda Moses 2018-12-14 Build site.
Rmd fd3d5ae Lambda Moses 2018-12-14 Don’t collapse output
html 2da3b35 Lambda Moses 2018-12-14 Build site.
Rmd 1fc3e91 Lambda Moses 2018-12-14 Publish 10xv2 notebook
Rmd 074d55f Lambda Moses 2018-12-14 Initial commit, already with 10xv2 notebook

In this vignette, we process fastq data from scRNA-seq (10x v2 chemistry) with to make a sparse matrix that can be used in downstream analysis. Then we will start that standard downstream analysis with Seurat.

Setup

First, if you want to run this notebook, please make sure you have git installed, go to the terminal and

git clone https://github.com/BUStools/BUS_notebooks_R.git

This will download a directory called BUS_notebooks_R in the directory you are in when running the git command. Inside the BUS_notebooks_R directory, this notebook is in the analysis directory. Then create an R project in the BUS_notebooks_R directory. Also, in RStudio, go to Tools -> Global Options -> R Markdown, and set “Evaluate chunks in directory” to “Project”. The reason why is that this website was built with workflowr, which renders the notebooks in the project directory rather than the document directory.

Alternatively, you can directly download this notebook here and open it in RStudio, but make sure to adjust the paths where files are downloaded and where output files are stored to match those in your environment.

Install packages

The primary analysis section of this notebook demonstrates the use of command line tools kallisto and bustools. Please use kallisto 0.45, 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 here is the directory of interest. Then you can directly invoke the binary on the command line as we will do in this notebook.

Note for Windows users: bustools does not have a Windows binary, but you can use a Linux subsystem in Windows 10. If you use earlier version of Windows, then unfortunately you need to either use Linux dual boot or virtual machine or find a Linux or Mac computer such as a server.

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. We will also use Seurat 3.0, which is not yet on CRAN (the CRAN version is 2.3.4), for the standard downstream analysis. Seurat 3.0 has some cool new features such as getting batch corrected gene count matrix and cell label transfer across datasets. See this paper for more detail. Though we will not get into those new features here, it’s worthwhile to check them out.

# 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", ref = "release/3.0")

The package DropletUtils will be used to estimate the number of real cells as opposed to empty droplets. It’s on Bioconductor, and here is how it should be installed:

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

The other R packages below are on CRAN, and can be installed with install.packages.

library(BUSpaRse)
library(data.table)
library(Seurat)
library(tidyverse)
library(DropletUtils)
library(Matrix)

Download data

The data set we are using here is 1k 1:1 Mixture of Fresh Frozen Human (HEK293T) and Mouse (NIH3T3) Cells from the 10x website. First, we download the fastq files (6.34 GB).

if (!file.exists("./data/hgmm1k_fastqs.tar")) {
  download.file("http://cf.10xgenomics.com/samples/cell-exp/2.1.0/hgmm_1k/hgmm_1k_fastqs.tar", destfile = "./data/hgmm_1k_fastqs.tar", quiet = TRUE)
}

Then untar this file

cd ./data
tar -xvf ./hgmm_1k_fastqs.tar
#> fastqs/
#> fastqs/hgmm_1k_S1_L001_I1_001.fastq.gz
#> fastqs/hgmm_1k_S1_L001_R1_001.fastq.gz
#> fastqs/hgmm_1k_S1_L001_R2_001.fastq.gz
#> fastqs/hgmm_1k_S1_L002_I1_001.fastq.gz
#> fastqs/hgmm_1k_S1_L002_R1_001.fastq.gz
#> fastqs/hgmm_1k_S1_L002_R2_001.fastq.gz
#> fastqs/hgmm_1k_S1_L003_I1_001.fastq.gz
#> fastqs/hgmm_1k_S1_L003_R1_001.fastq.gz
#> fastqs/hgmm_1k_S1_L003_R2_001.fastq.gz
#> fastqs/hgmm_1k_S1_L004_I1_001.fastq.gz
#> fastqs/hgmm_1k_S1_L004_R1_001.fastq.gz
#> fastqs/hgmm_1k_S1_L004_R2_001.fastq.gz
#> fastqs/hgmm_1k_S1_L005_I1_001.fastq.gz
#> fastqs/hgmm_1k_S1_L005_R1_001.fastq.gz
#> fastqs/hgmm_1k_S1_L005_R2_001.fastq.gz
#> fastqs/hgmm_1k_S1_L006_I1_001.fastq.gz
#> fastqs/hgmm_1k_S1_L006_R1_001.fastq.gz
#> fastqs/hgmm_1k_S1_L006_R2_001.fastq.gz
#> fastqs/hgmm_1k_S1_L007_I1_001.fastq.gz
#> fastqs/hgmm_1k_S1_L007_R1_001.fastq.gz
#> fastqs/hgmm_1k_S1_L007_R2_001.fastq.gz
#> fastqs/hgmm_1k_S1_L008_I1_001.fastq.gz
#> fastqs/hgmm_1k_S1_L008_R1_001.fastq.gz
#> fastqs/hgmm_1k_S1_L008_R2_001.fastq.gz

Primary analysis

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 transcriptome. This data set has both human and mouse cells, so we need both human and mouse transcriptomes. The transcriptomes downloaded here are from Ensembl version 95, released in January 2019. As of the writing of this notebook, this is the most recent Ensembl release.

# Human transcriptome
if (!file.exists("./data/hs_cdna.fa.gz")) {
  download.file("ftp://ftp.ensembl.org/pub/release-95/fasta/homo_sapiens/cdna/Homo_sapiens.GRCh38.cdna.all.fa.gz", "./data/hs_cdna.fa.gz", method = "wget", quiet = TRUE)
}
# Mouse transcriptome
if (!file.exists("./data/mm_cdna.fa.gz")) {
  download.file("ftp://ftp.ensembl.org/pub/release-95/fasta/mus_musculus/cdna/Mus_musculus.GRCm38.cdna.all.fa.gz", "./data/mm_cdna.fa.gz", method = "wget", quiet = TRUE)
}
kallisto version
#> kallisto, version 0.45.0

Actually, we don’t need to unzip the fasta files

kallisto index -i ./output/hs_mm_tr_index.idx ./data/hs_cdna.fa.gz ./data/mm_cdna.fa.gz
#> 
#> [build] loading fasta file ./data/hs_cdna.fa.gz
#> [build] loading fasta file ./data/mm_cdna.fa.gz
#> [build] k-mer length: 31
#> [build] warning: clipped off poly-A tail (longer than 10)
#>         from 2083 target sequences
#> [build] warning: replaced 8 non-ACGUT characters in the input sequence
#>         with pseudorandom nucleotides
#> [build] counting k-mers ... done.
#> [build] building target de Bruijn graph ...  done 
#> [build] creating equivalence classes ...  done
#> [build] target de Bruijn graph has 2145397 contigs and contains 206575877 k-mers

Run kallisto bus

Here we will generate the bus file. Here bus stands for Barbode, UMI, Set (i.e. equivalent class). In text form, it is a table whose first column is the barcode. The second column is the UMI that are associated with the barcode. The third column is the index of the equivalence class reads with the UMI maps to (equivalence class will be explained later). The fourth column is count of reads with this barcode, UMI, and equivalence class combination, which is ignored as one UMI should stand for one molecule. See this paper for more detail.

These are the technologies supported by kallisto bus:

system("kallisto bus --list", intern = TRUE)
#> Warning in system("kallisto bus --list", intern = TRUE): running command
#> 'kallisto bus --list' had status 1
#>  [1] "List of supported single cell technologies"
#>  [2] ""                                          
#>  [3] "short name       description"              
#>  [4] "----------       -----------"              
#>  [5] "10Xv1            10X chemistry version 1"  
#>  [6] "10Xv2            10X chemistry verison 2"  
#>  [7] "DropSeq          DropSeq"                  
#>  [8] "inDrop           inDrop"                   
#>  [9] "CELSeq           CEL-Seq"                  
#> [10] "CELSeq2          CEL-Seq version 2"        
#> [11] "SCRBSeq          SCRB-Seq"                 
#> [12] ""                                          
#> attr(,"status")
#> [1] 1

Here we have 8 samples. Each sample has 3 files: I1 means sample index, R1 means barcode and UMI, and R2 means the piece of cDNA. The -i argument specifies the index file we just built. The -o argument specifies the output directory. The -x argument specifies the sequencing technology used to generate this data set. The -t argument specifies the number of threads used. I ran this on a server and used 8 threads.

cd ./data
kallisto bus -i ../output/hs_mm_tr_index.idx -o ../output/out_hgmm1k -x 10xv2 -t8 \
./fastqs/hgmm_1k_S1_L001_R1_001.fastq.gz ./fastqs/hgmm_1k_S1_L001_R2_001.fastq.gz \
./fastqs/hgmm_1k_S1_L002_R1_001.fastq.gz ./fastqs/hgmm_1k_S1_L002_R2_001.fastq.gz \
./fastqs/hgmm_1k_S1_L003_R1_001.fastq.gz ./fastqs/hgmm_1k_S1_L003_R2_001.fastq.gz \
./fastqs/hgmm_1k_S1_L004_R1_001.fastq.gz ./fastqs/hgmm_1k_S1_L004_R2_001.fastq.gz \
./fastqs/hgmm_1k_S1_L005_R1_001.fastq.gz ./fastqs/hgmm_1k_S1_L005_R2_001.fastq.gz \
./fastqs/hgmm_1k_S1_L006_R1_001.fastq.gz ./fastqs/hgmm_1k_S1_L006_R2_001.fastq.gz \
./fastqs/hgmm_1k_S1_L007_R1_001.fastq.gz ./fastqs/hgmm_1k_S1_L007_R2_001.fastq.gz \
./fastqs/hgmm_1k_S1_L008_R1_001.fastq.gz ./fastqs/hgmm_1k_S1_L008_R2_001.fastq.gz
#> 
#> [index] k-mer length: 31
#> [index] number of targets: 303,693
#> [index] number of k-mers: 206,575,877
#> [index] number of equivalence classes: 1,256,267
#> [quant] will process sample 1: ./fastqs/hgmm_1k_S1_L001_R1_001.fastq.gz
#>                                ./fastqs/hgmm_1k_S1_L001_R2_001.fastq.gz
#> [quant] will process sample 2: ./fastqs/hgmm_1k_S1_L002_R1_001.fastq.gz
#>                                ./fastqs/hgmm_1k_S1_L002_R2_001.fastq.gz
#> [quant] will process sample 3: ./fastqs/hgmm_1k_S1_L003_R1_001.fastq.gz
#>                                ./fastqs/hgmm_1k_S1_L003_R2_001.fastq.gz
#> [quant] will process sample 4: ./fastqs/hgmm_1k_S1_L004_R1_001.fastq.gz
#>                                ./fastqs/hgmm_1k_S1_L004_R2_001.fastq.gz
#> [quant] will process sample 5: ./fastqs/hgmm_1k_S1_L005_R1_001.fastq.gz
#>                                ./fastqs/hgmm_1k_S1_L005_R2_001.fastq.gz
#> [quant] will process sample 6: ./fastqs/hgmm_1k_S1_L006_R1_001.fastq.gz
#>                                ./fastqs/hgmm_1k_S1_L006_R2_001.fastq.gz
#> [quant] will process sample 7: ./fastqs/hgmm_1k_S1_L007_R1_001.fastq.gz
#>                                ./fastqs/hgmm_1k_S1_L007_R2_001.fastq.gz
#> [quant] will process sample 8: ./fastqs/hgmm_1k_S1_L008_R1_001.fastq.gz
#>                                ./fastqs/hgmm_1k_S1_L008_R2_001.fastq.gz
#> [quant] finding pseudoalignments for the reads ... done
#> [quant] processed 63,252,296 reads, 52,231,469 reads pseudoaligned

See what the outputs are

list.files("./output/out_hgmm1k/")
#> [1] "matrix.ec"         "output.bus"        "output.sorted"    
#> [4] "output.sorted.txt" "run_info.json"     "transcripts.txt"

Run BUStools

The output.bus file is a binary. In order to make R parse it, we need to convert it into a sorted text file. There’s a command line tool bustools for this. The reason why kallisto bus outputs this binary rather than the sorted text file directly is that the binary can be very compressed to save disk space. The such compression is achieved, a future version of R package BUSpaRse may support directly reading the binary to convert to sparse matrix. At present, as bustools is still under development, such compression has not been achieved, and as a result, the binary is slightly larger than the text file.

# Add where I installed bustools to PATH
export PATH=$PATH:/home/lambda/mylibs/bin/
# Sort, with 8 threads
bustools sort -o ./output/out_hgmm1k/output.sorted -t8 ./output/out_hgmm1k/output.bus
# Convert sorted file to text
bustools text -o ./output/out_hgmm1k/output.sorted.txt ./output/out_hgmm1k/output.sorted
#> Read in 52231469 number of busrecords
#> All sorted
#> Read in 43778790 number of busrecords

This is what the bus file looks like:

First column is barcode, second column is UMI, 3rd column is equivalence class index, and 4th column is count, which is ignored. This has already been sorted according to barcode.

head ./output/out_hgmm1k/output.sorted.txt
#> AAAAAAAAAAAAAAAA AAAAAAAAAA  6575    1
#> AAAAAAAAAAAAAAAA AAAAAAAAAA  9414    1
#> AAAAAAAAAAAAAAAA AAAAAAAAAA  10354   1
#> AAAAAAAAAAAAAAAA AAAAAAAAAA  38532   1
#> AAAAAAAAAAAAAAAA AAAAAAAAAA  112639  1
#> AAAAAAAAAAAAAAAA AAAAAAAAAA  123866  1
#> AAAAAAAAAAAAAAAA AAAAAAAAAA  136025  1
#> AAAAAAAAAAAAAAAA AAAAAAAAAA  204499  2
#> AAAAAAAAAAAAAAAA AAAAAAAAAA  209310  1
#> AAAAAAAAAAAAAAAA AAAAAAAAAA  221957  1

Secondary analysis

Map transcripts to genes

For the sparse matrix, most people are interested in how many UMIs per gene per cell, we here we will quantify this from the bus output, and to do so, we need to find which gene corresponds to each transcript. Remember in the output of kallisto bus, there’s the file transcripts.txt. Those are the transcripts in the transcriptome index. Now we’ll only keep the transcripts present there and make sure that the transcripts in tr2g are in the same order as those in the index. You will soon see why the order is important.

Remember that we downloaded transcriptome FASTA files from Ensembl just now. In FASTA files, each entry is a sequence with a name. In Ensembl FASTA files, the sequence name has genome annotation of the corresponding sequence, so we can extract transcript IDs and corresponding gene IDs and gene names from there.

tr2g <- transcript2gene(fasta_file = c("./data/hs_cdna.fa.gz", "./data/mm_cdna.fa.gz"),
                        kallisto_out_path = "./output/out_hgmm1k")
#> Reading FASTA file.
#> Reading FASTA file.
#> Sorting transcripts
head(tr2g) %>% knitr::kable()
transcript gene gene_name
ENST00000632684.1 ENSG00000282431.1 TRBD1
ENST00000434970.2 ENSG00000237235.2 TRDD2
ENST00000448914.1 ENSG00000228985.1 TRDD3
ENST00000415118.1 ENSG00000223997.1 TRDD1
ENST00000631435.1 ENSG00000282253.1 TRBD1
ENST00000390583.1 ENSG00000211923.1 IGHD3-10

Alternative ways of getting tr2g have been implemented in the BUSpaRse package. You may use tr2g_ensembl to query Ensembl with biomart to get transcript and gene IDs. If you use this method, then please make sure that the Ensembl version used in the query matches that of the transcriptome. This method is convenient for the user since you only need to input species names, but it can be slow since biomart database query can be slow. You may also use tr2g_gtf for GTF files and tr2g_gff3 for GFF3 files, which are more useful for non-model organisms absent from Ensemble. After calling the tr2g_* family of functions, you should sort the transcripts from those functions with sort_tr2g so the transcripts are in the same order as those in the kallisto index. The FASTA way is the fastest, as reading FASTA files is faster than reading GTF and GFF files. Whichever way you choose, this shouldn’t take much more than half a minute per species.

The function transcript2gene not only gets the transcript and gene IDs but also sorts the transcripts. This function only supports Ensembl biomart query and Ensembl fasta files as the source of transcript and gene IDs, as the attribute field of GTF and GFF3 files can differ between sources and further cleanup may be needed.

Make the sparse matrix

For 10x, we do have a file with all valid cell barcodes that comes with CellRanger. You need to install CellRanger to get this file, though you do not need to run CellRanger for this notebook. The whitelist is optional, so if you don’t have one, you may skip the whitelist step and the whitelist argument in the makr_sparse_matrix function.

# Copy v2 chemistry whitelist to working directory
cp ~/cellranger-3.0.1/cellranger-cs/3.0.1/lib/python/cellranger/barcodes/737K-august-2016.txt \
./data/whitelist_v2.txt
# Read in the whitelist
whitelist_v2 <- fread("./data/whitelist_v2.txt", header = FALSE)$V1
length(whitelist_v2)
#> [1] 737280

There about 737K valid cell barcodes.

Now we have everything we need to make the sparse matrix. This function reads in output.sorted.txt line by line and processes them. It does not do barcode correction for now, so the barcode must exactly match those in the whitelist if one is provided. It took 5 to 6 minutes to construct the sparse matrix in the hgmm6k dataset, which has over 280 million lines in output.sorted.txt, which is over 9GB. Here the data set is smaller, and it takes less than a minute. A future version of BUSpaRse will parallelize this process to make it more scalable for huge datasets. For now, it’s single threaded.

Note that the arguments est_ncells (estimated number of cells) and est_ngenes (estimated number of genes) are important. With the estimate, this function reserves memory for the data to be added into, reducing the need of reallocation, which will slow the function down. Since the vast majority of “cells” you get in this sparse matrix are empty droplets rather than cells, please put at least 200 times more “cells” than you actually expect in est_ncells.

If you do not have a whitelist of barcodes, then it’s fine; the whitelist argument is optional.

The function make_sparse_matrix can make the gene count matrix and the transcript compatibility count (TCC) matrix at the same time. For the purpose of this notebook, we only generate the gene count matrix. An upcoming notebook will demonstrate some more detailed analysis with a TCC matrix. See Ntranos et al. 2016 for more information about TCC matrices.

res_mat <- make_sparse_matrix("./output/out_hgmm1k/output.sorted.txt",
                              tr2g = tr2g, est_ncells = 3.5e5,
                              est_ngenes = nrow(tr2g),
                              whitelist = whitelist_v2, TCC = FALSE)
#> Reading matrix.ec
#> Processing ECs
#> Matching genes to ECs
#> Reading data
#> Read 5 million reads
#> Read 10 million reads
#> Read 15 million reads
#> Read 20 million reads
#> Read 25 million reads
#> Read 30 million reads
#> Read 35 million reads
#> Read 40 million reads
#> Constructing gene count matrix

The matrix we get has genes in rows and barcode in columns. The row names are the gene IDs (not using human readable gene names since they’re not guaranteed to be unique), and the column names are cell barcodes.

For Ensembl transcriptomes, you may call busparse_matrix to bypass the separate step of getting the tr2g data frame here to directly get the matrix.

Explore the data

Remove empty droplets

Cool, so now we have the sparse matrix. What does it look like?

dim(res_mat)
#> [1]  47367 353961

That’s way more cells than we expect, which is about 1000. So what’s going on?

How many UMIs per barcode?

tot_counts <- Matrix::colSums(res_mat)
summary(tot_counts)
#>    Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
#>     1.0     1.0     2.0   101.8    11.0 92917.0

The vast majority of “cells” have only 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 method 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)
qplot(bc_rank$rank, bc_rank$total, geom = "line") +
  geom_hline(yintercept = bc_rank$knee, color = "blue", linetype = 2) +
  geom_hline(yintercept = bc_rank$inflection, color = "green", linetype = 2) +
  annotate("text", x = 1000, y = 1.5 * c(bc_rank$knee, bc_rank$inflection),
           label = c("knee", "inflection"), color = c("blue", "green")) +
  scale_x_log10() +
  scale_y_log10() +
  labs(x = "Rank", y = "Total reads") +
  theme_bw()

Version Author Date
8cd81e0 Lambda Moses 2019-02-14

The inflection point looks like a reasonable number of cells.

# Filter the matrix
res_mat <- res_mat[, tot_counts > bc_rank$inflection]
dim(res_mat)
#> [1] 47367  1123

Cell species

How many cells are from humans and how many from mice? The number of cells with mixed species indicates doublet rate.

gene_species <- ifelse(str_detect(rownames(res_mat), "^ENSMUSG"), "mouse", "human")
mouse_inds <- gene_species == "mouse"
human_inds <- gene_species == "human"
# mark cells as mouse or human
cell_species <- tibble(n_mouse_umi = Matrix::colSums(res_mat[mouse_inds,]),
                       n_human_umi = Matrix::colSums(res_mat[human_inds,]),
                       tot_umi = Matrix::colSums(res_mat),
                       prop_mouse = n_mouse_umi / tot_umi,
                       prop_human = n_human_umi / tot_umi)
# Classify species based on proportion of UMI, with cutoff of 90%
cell_species <- cell_species %>% 
  mutate(species = case_when(
    prop_mouse > 0.9 ~ "mouse",
    prop_human > 0.9 ~ "human",
    TRUE ~ "mixed"
  ))
ggplot(cell_species, aes(n_human_umi, n_mouse_umi, color = species)) +
  geom_point(size = 0.5) +
  theme_bw()

Version Author Date
8cd81e0 Lambda Moses 2019-02-14

Great, looks like the vast majority of cells are not mixed.

cell_species %>% 
  dplyr::count(species) %>% 
  mutate(proportion = n / ncol(res_mat))
#> # A tibble: 3 x 3
#>   species     n proportion
#>   <chr>   <int>      <dbl>
#> 1 human     579    0.516  
#> 2 mixed       3    0.00267
#> 3 mouse     541    0.482

Great, only about 0.3% of cells here are doublets, which is lower than the ~1% 10x lists. Doublet rate tends to be lower when cell concentration is lower. However, doublets can still be formed with cells from the same species, so the number of mixed species “cells” is only a lower bound of doublet rate.

Dimension reduction

seu <- CreateSeuratObject(res_mat, min.cells = 3) %>% 
  NormalizeData(verbose = FALSE) %>% 
  ScaleData(verbose = FALSE) %>% 
  FindVariableFeatures(verbose = FALSE)
# Add species to meta data
seu <- AddMetaData(seu, metadata = cell_species$species, col.name = "species")

See how number of total counts and number of genes expressed are distributed.

VlnPlot(seu, c("nCount_RNA", "nFeature_RNA"), group.by = "species")

Version Author Date
8cd81e0 Lambda Moses 2019-02-14
seu <- RunPCA(seu, verbose = FALSE, npcs = 30)
ElbowPlot(seu, ndims = 30)

Version Author Date
8cd81e0 Lambda Moses 2019-02-14
DimPlot(seu, reduction = "pca", pt.size = 0.5, group.by = "species")

Version Author Date
8cd81e0 Lambda Moses 2019-02-14

The first PC separates species, as expected. Also as expected, the doublets are in between human and mouse cells in this plot.

seu <- RunTSNE(seu, dims = 1:20, check_duplicates = FALSE)
DimPlot(seu, reduction = "tsne", pt.size = 0.5, group.by = "species")

Version Author Date
8cd81e0 Lambda Moses 2019-02-14
The species separate, and the few doublets form its own cluster, as expected.

sessionInfo()
#> R version 3.5.1 (2018-07-02)
#> Platform: x86_64-redhat-linux-gnu (64-bit)
#> Running under: CentOS Linux 7 (Core)
#> 
#> Matrix products: default
#> BLAS/LAPACK: /usr/lib64/R/lib/libRblas.so
#> 
#> locale:
#>  [1] LC_CTYPE=en_US.UTF-8       LC_NUMERIC=C               LC_TIME=en_US.UTF-8       
#>  [4] LC_COLLATE=en_US.UTF-8     LC_MONETARY=en_US.UTF-8    LC_MESSAGES=en_US.UTF-8   
#>  [7] LC_PAPER=en_US.UTF-8       LC_NAME=C                  LC_ADDRESS=C              
#> [10] LC_TELEPHONE=C             LC_MEASUREMENT=en_US.UTF-8 LC_IDENTIFICATION=C       
#> 
#> attached base packages:
#> [1] parallel  stats4    stats     graphics  grDevices utils     datasets  methods  
#> [9] base     
#> 
#> other attached packages:
#>  [1] Matrix_1.2-15               DropletUtils_1.2.2          SingleCellExperiment_1.4.1 
#>  [4] SummarizedExperiment_1.12.0 DelayedArray_0.8.0          matrixStats_0.54.0         
#>  [7] Biobase_2.42.0              GenomicRanges_1.34.0        GenomeInfoDb_1.18.1        
#> [10] IRanges_2.16.0              S4Vectors_0.20.1            BiocGenerics_0.28.0        
#> [13] BiocParallel_1.16.5         forcats_0.4.0               stringr_1.4.0              
#> [16] dplyr_0.8.0.1               purrr_0.3.0                 readr_1.3.1                
#> [19] tidyr_0.8.2                 tibble_2.0.1                ggplot2_3.1.0              
#> [22] tidyverse_1.2.1             Seurat_3.0.0.9000           data.table_1.12.0          
#> [25] BUSpaRse_0.99.0             workflowr_1.2.0            
#> 
#> loaded via a namespace (and not attached):
#>   [1] readxl_1.3.0             backports_1.1.3          plyr_1.8.4              
#>   [4] igraph_1.2.4             lazyeval_0.2.1           splines_3.5.1           
#>   [7] listenv_0.7.0            usethis_1.4.0            digest_0.6.18           
#>  [10] htmltools_0.3.6          fansi_0.4.0              gdata_2.18.0            
#>  [13] magrittr_1.5             memoise_1.1.0            cluster_2.0.7-1         
#>  [16] ROCR_1.0-7               limma_3.38.3             remotes_2.0.2           
#>  [19] globals_0.12.4           Biostrings_2.50.2        modelr_0.1.4            
#>  [22] RcppParallel_4.4.2       R.utils_2.8.0            prettyunits_1.0.2       
#>  [25] colorspace_1.4-0         rvest_0.3.2              blob_1.1.1              
#>  [28] ggrepel_0.8.0            haven_2.1.0              xfun_0.5                
#>  [31] callr_3.1.1              crayon_1.3.4             RCurl_1.95-4.11         
#>  [34] jsonlite_1.6             zeallot_0.1.0            survival_2.43-3         
#>  [37] zoo_1.8-4                ape_5.2                  glue_1.3.0              
#>  [40] gtable_0.2.0             zlibbioc_1.28.0          XVector_0.22.0          
#>  [43] pkgbuild_1.0.2           plyranges_1.2.0          Rhdf5lib_1.4.0          
#>  [46] future.apply_1.1.0       HDF5Array_1.10.1         scales_1.0.0            
#>  [49] edgeR_3.24.3             DBI_1.0.0                bibtex_0.4.2            
#>  [52] Rcpp_1.0.0               metap_1.1                viridisLite_0.3.0       
#>  [55] progress_1.2.0           reticulate_1.11          bit_1.1-14              
#>  [58] rsvd_1.0.0               SDMTools_1.1-221         tsne_0.1-3              
#>  [61] htmlwidgets_1.3          httr_1.4.0               gplots_3.0.1.1          
#>  [64] RColorBrewer_1.1-2       ica_1.0-2                pkgconfig_2.0.2         
#>  [67] XML_3.98-1.17            R.methodsS3_1.7.1        utf8_1.1.4              
#>  [70] locfit_1.5-9.1           labeling_0.3             tidyselect_0.2.5        
#>  [73] rlang_0.3.1              AnnotationDbi_1.44.0     cellranger_1.1.0        
#>  [76] munsell_0.5.0            tools_3.5.1              cli_1.0.1               
#>  [79] generics_0.0.2           RSQLite_2.1.1            broom_0.5.1             
#>  [82] devtools_2.0.1           ggridges_0.5.1           evaluate_0.13           
#>  [85] yaml_2.2.0               npsurv_0.4-0             processx_3.2.1          
#>  [88] knitr_1.21               bit64_0.9-7              fs_1.2.6                
#>  [91] fitdistrplus_1.0-14      caTools_1.17.1.1         RANN_2.6.1              
#>  [94] pbapply_1.4-0            future_1.11.1.1          nlme_3.1-137            
#>  [97] whisker_0.3-2            R.oo_1.22.0              xml2_1.2.0              
#> [100] biomaRt_2.38.0           rstudioapi_0.9.0         compiler_3.5.1          
#> [103] plotly_4.8.0             png_0.1-7                testthat_2.0.1          
#> [106] lsei_1.2-0               stringi_1.3.1            highr_0.7               
#> [109] ps_1.3.0                 desc_1.2.0               lattice_0.20-38         
#> [112] pillar_1.3.1             Rdpack_0.10-1            lmtest_0.9-36           
#> [115] cowplot_0.9.4            bitops_1.0-6             irlba_2.3.3             
#> [118] gbRd_0.4-11              rtracklayer_1.42.1       R6_2.4.0                
#> [121] KernSmooth_2.23-15       sessioninfo_1.1.1        codetools_0.2-16        
#> [124] MASS_7.3-51.1            gtools_3.8.1             assertthat_0.2.0        
#> [127] pkgload_1.0.2            rhdf5_2.26.0             rprojroot_1.3-2         
#> [130] withr_2.1.2              GenomicAlignments_1.18.1 Rsamtools_1.34.1        
#> [133] GenomeInfoDbData_1.2.0   hms_0.4.2                grid_3.5.1              
#> [136] rmarkdown_1.11           Rtsne_0.15               git2r_0.24.0            
#> [139] lubridate_1.7.4