你想在笔记本上跑CellRanger么?
分享是一种态度
这个题目有点做梦。但是疫情这么严重，我从GEO下载的单细胞RNA-seq(scRNA-seq)原始数据想做下游分析，手头只有一个2015年已经热变形的MacBookPro （2.5GHz i7 的CPU 16G RAM的内存），我实在不忍心再伤害它了。除了CellRanger，还有没有什么替代办法能在笔记本上跑呢？我想到了自己2016年就开始用的Kallisto，因为能耗小，速度快，我经常在自己笔记本上跑RNA-seq，一个样平均半小时的速度让我从入门RNA-seq就放弃了Bowtie2+TopHat或者STAR。那么，这个程序能不能跑scRNA-seq呢?我翻到了前段时间的溜得飞起的单细胞bus你还不上车？,顺着题目找到了官网https://www.kallistobus.tools/
这个主页的封面让我想起了kallisto（https://pachterlab.github.io/kallisto/）的主页图
意思很明确了，一群kallisto上车了~
简单来说Kallisto Bustools (后文简称KB)只有两套命令，ref和count，前者是利用cDNA序列和gtf注释文件建立index，而count则是输出表达矩阵。值得注意的是，ref命令输出内含子和剪切信息，因此KB可以生成适合RNA Velocyto的的矩阵信息。网站除了两个命令的说明外还有非常详细的实例说明，在科学上网的情况下可以直接在colab环境测试运行。
1 KB的安装
1.1 环境准备
硬件方面，如开头所述，我的电脑是2015年中的MBP15‘，非常普通而且老旧的笔记本电脑。操作系统是Catalina(10.15.3)，Xcode （11.3.1）。
首先使用Anaconda建立虚拟环境，这里有一点我栽过跟头，建立环境时一定注明Python版本，否则安装kb时会有冲突。
conda create --name KB_kite python=3.8 #我的系统有俩python3一个3.7一个3.8，如果只写3，后面安装会报错。
conda activate KB_kite
这里注意我的环境名有意写成KB_kite而不是KB是因为我实际上安装了两个版本的KB，如果后面需要使用RNA Velocyto，建议安装count_kite，否则无法建立正确的index。
1.2 安装
KB的安装就是从github获得程序进行安装，而1.1所提到的count_kite是KB的一个branch。安装有两种方式：
方式一：
pip install git+https://github.com/pachterlab/kb_python@count-kite #建议科学上网，否则很容易报错
方式二：
直接登录https://github.com/pachterlab/kb_python/tree/count-kite下载，然后进行本地安装：
pip install ./kb_python-count-kite.zip
安装需要大概十来分钟就结束。
2 索引建立
如前所述，索引（index）建立只用一个ref就完成。有很多的参数，具体可以参考官网，值得注意的是如果是要建立RNA Velocyto相关的索引，那么一定要使用count-kite分支，否则lamanno不会被识别。此外，还有个选项-d是直接下载官网建立好的索引数据，如果网络不好，慎用！
这里我根据我自己的工作选择了小鼠GRCm38的cDNA序列：
wget ftp://ftp.ensembl.org/pub/release-98/fasta/mus_musculus/dna/Mus_musculus.GRCm38.dna.primary_assembly.fa.gz
wget ftp://ftp.ensembl.org/pub/release-98/gtf/mus_musculus/Mus_musculus.GRCm38.98.gtf.gz
下载完后直接建立索引，根据官方说明，我需要为后续RNA Velocyto做准备，因此选择了--workflow lamanno，具体命令如下：
kb ref -i index.idx -g t2g.txt -f1 cdna.fa -f2 intron.fa -c1 cdna_t2c.txt -c2 intron_t2c.txt --workflow lamanno -n 8 \
Mus_musculus.GRCm38.dna.primary_assembly.fa.gz \
Mus_musculus.GRCm38.98.gtf.gz #-n 8是将最终的index分成8份，这样减少了我电脑的压力
我记得以前用STAR给人cDNA做索引时32G内存折腾了快一天，电脑假死了很久，这次的索引建立，请从下面截图感受一下：
电脑的风扇哀嚎了一小时就完成。最终产出的文件有8个index.idx文件，cdna.fa，cdna_t2c.txt，intron.fa，intron_t2c.txt和基因名注释文件：t2g.txt
3 表达矩阵产生
3.1 测序数据准备
这里为了广大CellRanger的学员比较时间，我采用的是10x官方数据pbmc_1k，并用上述命令建立了人转录组的索引。
wget http://cf.10xgenomics.com/samples/cell-exp/3.0.0/pbmc_1k_v3/pbmc_1k_v3_fastqs.tar #下载数据
tar -xvf pbmc_1k_v3_fastqs.tar
3.2 表达矩阵产生
非常简单的命令，注意参数-x 根据建库试剂盒填写，这里的数据是10xv3。
kb count -i index.idx -g t2g.txt -x 10xv3 -o output --filter bustools -t 2 pbmc_1k_v3_fastqs/pbmc_1k_v3_S1_L001_R1_001.fastq.gz pbmc_1k_v3_fastqs/pbmc_1k_v3_S1_L001_R2_001.fastq.gz pbmc_1k_v3_fastqs/pbmc_1k_v3_S1_L002_R1_001.fastq.gz pbmc_1k_v3_fastqs/pbmc_1k_v3_S1_L002_R2_001.fastq.gz
结果会在output文件夹产生两个子文件夹，一个是filtered一个是unfiltered，每个子文件夹就三个数据：cells_x_genes.barcodes.txt，cells_x_genes.genes.txt，cells_x_genes.mtx，熟悉10x数据的人看名字就看得懂这三个各代表啥了。
哦，我忘记截图了，整个表达矩阵产生时间是17分23秒。
3.3 常用的QC
这里选择unfiltered的三个文件进行质检，这里参考了官方的R笔记（https://colab.research.google.com/github/pachterlab/kallistobustools/blob/master/notebooks/kb_analysis_0_R.ipynb#scrollTo=Z3MgRWtoTbL0） (我知道你们一般点了进不去，所以我一边搬运一边学习了)
R包的调取
library(DropletUtils)
library(Seurat)
library(Matrix)
library(tidyverse)
library(scico)
setwd(“path/to/your/output”)
读取表达矩阵(我觉得这个函数很精巧的)
read_count_output <- function(dir, name) {
dir <- normalizePath(dir, mustWork = TRUE)
m <- readMM(paste0(dir, "/", name, ".mtx"))
m <- Matrix::t(m)
m <- as(m, "dgCMatrix")
# The matrix read has cells in rows
ge <- ".genes.txt"
genes <- readLines(file(paste0(dir, "/", name, ge)))
barcodes <- readLines(file(paste0(dir, "/", name, ".barcodes.txt")))
colnames(m) <- barcodes
rownames(m) <- genes
return(m)
}
pbmc_mat <- read_count_output("./output/counts_unfiltered", name = "cells_x_genes") #读取表达矩阵
dim(pbmc_mat) #60623,259615
检查测序饱和度：
total_counts <- colSums(pbmc_mat)
pbmc_sat <- tibble(nCount = total_counts,
nGene = colSums(pbmc_mat > 0))
options(repr.plot.width=9, repr.plot.height=6)
ggplot(pbmc_sat, aes(nCount, nGene)) +
geom_bin2d(bins = 50) +
scale_fill_scico(palette = "devon", direction = -1, end = 0.95) +
scale_x_log10() + scale_y_log10() + annotation_logticks()
knee plot检查（umm…我也觉得这个函数很精巧所以拿来了）
knee_plot <- function(bc_rank) {
knee_plt <- tibble(rank = bc_rank[["rank"]],
total = bc_rank[["total"]]) %>%
distinct() %>%
dplyr::filter(total > 0)
annot <- tibble(inflection = metadata(bc_rank)[["inflection"]],
rank_cutoff = max(bc_rank$rank[bc_rank$total > metadata(bc_rank)[["inflection"]]]))
p <- ggplot(knee_plt, aes(total, rank)) +
geom_line() +
geom_hline(aes(yintercept = rank_cutoff), data = annot, linetype = 2) +
geom_vline(aes(xintercept = inflection), data = annot, linetype = 2) +
scale_x_log10() +
scale_y_log10() +
annotation_logticks() +
labs(y = "Rank", x = "Total UMIs")
return(p)
}
bc_rank <- barcodeRanks(pbmc_mat, lower = 1000)
options(repr.plot.width=9, repr.plot.height=6)
knee_plot(bc_rank)
根据上述结果进行简单的滤过空的droplet数据：
pbmc_mat <- pbmc_mat[, total_counts > metadata(bc_rank)$inflection]
pbmc_mat <- pbmc_mat[Matrix::rowSums(pbmc_mat) > 0,]
dim(res_mat) # 31832,1178
ensembleID向Gene Symbol转换：
t2g <- read_tsv("t2g.txt", col_names = c("transcript", "gene", "gene_name")) #t2g.txt为KB index产生文件
t2g <- distinct(t2g[, c("gene", "gene_name")])
rownames(pbmc_mat) <- t2g$gene_name[match(rownames(pbmc_mat), t2g$gene)]
3.4 数据导入Seurat
经过上述过滤（当然你还可以进行其他质检），我们导入Seurat
pbmc <- CreateSeuratObject(counts = pbmc_mat, project = "pbmc1k", min.cells = 3, min.features = 200) #后面的你都懂了吧
4 个人感受
我觉得KB保持了kallisto的一贯作风--“快，准“，没有CellRanger报告那么精美，但是给了非常实在的数字，QC方面可以自己用Rmarkdown做的更漂亮。虽然我没有在我自己电脑上尝试CellRanger，但是根据我以往STAR的使用经历，我不敢说KB很快，但是其消耗和速度对我和我的笔记本来说是非常友善。
如果你对单细胞转录组研究感兴趣，但又不知道如何入门，也许你可以关注一下下面的课程
看完记得顺手点个“在看”哦！
长按扫码可关注