<jupyter_start><jupyter_text><jupyter_code>install.packages("BiocManager")
BiocManager::install(c("multtest", "DropletUtils", "BUSpaRse", "TENxBUSData"))
install.packages(c("Seurat", "biomartr", "data.table"))
library(BUSpaRse)
library(TENxBUSData)
library(ggplot2)
library(magrittr)
library(data.table)
library(Seurat)
library(DropletUtils)
library(Matrix)
library(tibble)
theme_set(theme_bw())<jupyter_output><empty_output><jupyter_text>We will not get into the details of how to make the `bus` file with `kallisto bus` and `bustools`, as the data will be downloaded with `TENxBUSData`. `TENxBUSData` provides 5 different datasets, and we will use the PBMC 1k dataset here. The data from `TENxBUSData` contains the sorted `bus` file in text format. While the `BUSpaRse` package converts that text format `bus` file into gene count matrix, this text file can be loaded into R as a data frame for further exploration.<jupyter_code>fn <- TENxBUSData(".", dataset = "pbmc1k")
list.files(fn)<jupyter_output><empty_output><jupyter_text>Explaining the output:* `matrix.ec`: A text file with two columns. The first column is the 0 based index of equivalence classes. The second column is the set of transcripts (denoted by 0 based index based on order of appearance in the transcriptome fasta file) present in the corresponding equivalence class.* `output.sorted`: The data represented in `bus` format, sorted by barcode, UMI, and equivalence class. This is a binary, so can't be read into R with functions like `read.table`.* `output.sorted.txt`: `output.sorted` converted into text format, so can be easily read into R for exploration.* `transcript.txt`: A text file with one column, which is the transcripts present in the data, in the same order as in the transcriptome fasta file. Sparse matrix Map transcripts to genesFor 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. Information on which transcript corresponds to which gene can be directly retrieved from Ensembl.<jupyter_code>tr2g <- transcript2gene(species = "Homo sapiens", 
                        kallisto_out_path = "./out_pbmc1k",
                        ensembl_version = 94)
head(tr2g)<jupyter_output><empty_output><jupyter_text>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. Then the function `save_tr2g_bustools` can be used to save the `tr2g` data frame to a text file in the format required by `bustools`. Make the sparse matrixFor 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.<jupyter_code>if (!dir.exists("./reference")) dir.create("./reference")
download.file("https://caltech.box.com/public/static/q078m83pq1yumorcsbsom7en6m6uonqa.gz", 
destfile = "./reference/whitelist_v3.txt.gz", method = "wget")
# Read in the whitelist
whitelist_v3 <- fread("./reference/whitelist_v3.txt.gz", header = FALSE)$V1
length(whitelist_v3)<jupyter_output><empty_output><jupyter_text>That's an order of magnitude more than the 737K in v2 chemistry.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. 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](https://genomebiology.biomedcentral.com/articles/10.1186/s13059-016-0970-8) for more information about TCC matrices.<jupyter_code>res_mat <- make_sparse_matrix("./out_pbmc1k/output.sorted.txt",
                              tr2g = tr2g, est_ncells = 3e5,
                              est_ngenes = nrow(tr2g),
                              whitelist = whitelist_v3, TCC = FALSE)<jupyter_output>Reading matrix.ec
Processing ECs
Matching genes to ECs
Reading data
Read 5 million reads
Read 10 million reads
Read 15 million reads
Constructing gene count matrix<jupyter_text>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. Explore the data Remove empty dropletsCool, so now we have the sparse matrix. What does it look like?<jupyter_code>dim(res_mat)<jupyter_output><empty_output><jupyter_text>That's way more cells than we expect, which is about 1000. So what's going on?How many UMIs per barcode?<jupyter_code>tot_counts <- Matrix::colSums(res_mat)
summary(tot_counts)<jupyter_output><empty_output><jupyter_text>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`](https://www.bioconductor.org/packages/devel/bioc/vignettes/DropletUtils/inst/doc/DropletUtils.htmldetecting-empty-droplets)), for simplicity, we will use the barcode ranking method here. However, whichever way we go, we don't have the ground truth.<jupyter_code># Compute barcode rank
bc_rank <- barcodeRanks(res_mat)
#' Knee plot for filtering empty droplets
#' 
#' Visualizes the inflection point to filter empty droplets. This function plots 
#' different datasets with a different color. Facets can be added after calling
#' this function with `facet_*` functions.
#' 
#' @param bc_rank A `DataFrame` output from `DropletUtil::barcodeRanks`.
#' @return A ggplot2 object.
knee_plot <- function(bc_rank) {
  knee_plt <- data.frame(rank = bc_rank[["rank"]],
                         total = bc_rank[["total"]]) %>% 
    distinct() %>% 
    dplyr::filter(total > 0)
  annot <- data.frame(inflection = metadata(bc_rank)[["inflection"]],
                      rank_cutoff = max(bc_rank$rank[bc_rank$total > metadata(bc_rank)[["inflection"]]]))
  p <- ggplot(knee_plt, aes(rank, total)) +
    geom_line() +
    geom_hline(aes(yintercept = inflection), data = annot, linetype = 2) +
    geom_vline(aes(xintercept = rank_cutoff), data = annot, linetype = 2) +
    scale_x_log10() +
    scale_y_log10() +
    annotation_logticks() +
    labs(x = "Rank", y = "Total UMIs")
  return(p)
}
options(repr.plot.width=9, repr.plot.height=6)
knee_plot(bc_rank)<jupyter_output><empty_output><jupyter_text>The inflection point looks like a reasonable number of cells.<jupyter_code># Filter the matrix
res_mat <- res_mat[, tot_counts > metadata(bc_rank)$inflection]
res_mat <- res_mat[Matrix::rowSums(res_mat) > 0,]
dim(res_mat)<jupyter_output><empty_output><jupyter_text>Dimension reduction<jupyter_code>seu <- CreateSeuratObject(res_mat, min.cells = 3) %>% 
  SCTransform(verbose = FALSE)
VlnPlot(seu, c("nCount_RNA", "nFeature_RNA"))
ggplot(seu@<EMAIL>, aes(nCount_RNA, nFeature_RNA)) +
  geom_point() +
  labs(x = "Total UMI counts per cell", y = "Number of genes detected")
seu <- RunPCA(seu, verbose = FALSE, npcs = 30)
ElbowPlot(seu, ndims = 30)
seu <- FindNeighbors(seu)
seu <- FindClusters(seu)
PCAPlot(seu)
seu <- RunUMAP(seu, dims = 1:20, verbose = FALSE)
UMAPPlot(seu)<jupyter_output><empty_output>
{'fraction_non_alphanumeric': 0.08415426535317544, 'fraction_numerical': 0.011947180884510585, 'mean_word_length': 4.836697247706422, 'pattern_counts': {'":': 0, '<': 58, '<?xml version=': 0, '>': 48, 'https://': 3, 'lorem ipsum': 0, 'www.': 1, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 2, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
{'chain_length': '10', 'dir': 'jupyter-structured-clean-dedup', 'id': '614469', 'license': 'BSD-3-Clause', 'n_tokens_mistral': 3057, 'n_tokens_neox': 2817, 'n_words': 1155, 'path': 'notebooks/kb_10xv3.ipynb', 'repo_name': 'tar-singh/kallistobustools', 'size': '9544'}