你想在笔记本上跑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很快，但是其消耗和速度对我和我的笔记本来说是非常友善。