Kallisto Bustools

kallisto | bustools: Pseudoalignment-Based scRNA-seq Quantification

Overview

kallisto + bustools is a pseudoalignment-based quantification pipeline for single-cell RNA-seq. The standard tools (cellranger, STARsolo, alevin-fry) do full alignment of reads to the genome; kallisto skips that and only determines which transcripts are compatible with each read. The result: 5–20× faster than cellranger, comparable accuracy, much lower memory footprint.

The toolkit has three components:

For most users, kb-python is the only entry point you need. This protocol uses kb throughout.

When to Use This Skill

• Generating count matrices from raw FASTQ files for any droplet-based chemistry (10X, Drop-seq, inDrops) or plate-based (SMART-Seq)
• Processing many samples where cellranger would take days — kb does it in hours
• snRNA-seq / RNA velocity workflows that need both nascent (pre-mRNA) and mature counts
• Custom seqspec-defined chemistries
• Pipelines where you want maximum reproducibility — kb is deterministic, while cellranger has subtle non-determinism

Not for:

• Cell calling / empty-drop detection — kb produces unfiltered matrices; use DropletUtils::emptyDrops or CellBender downstream
• Cell-type annotation, clustering, DE — handoff to scanpy / Seurat / SCE
• Spatial transcriptomics — use Space Ranger or the spatial-transcriptomics protocol

Prerequisites

• Python 3.8+
• ~32 GB RAM for human/mouse index build; ~16 GB for kb count per sample
• Linux/macOS preferred (Windows works via WSL)
• Reference: a transcriptome FASTA + GTF (Ensembl, GENCODE, RefSeq)

Installation

# pip (recommended)
pip install kb_python

# OR conda
conda install -c bioconda kb-python

# Verify
kb --version
kb --list                              # supported single-cell technologies

kb-python bundles the kallisto + bustools binaries — no separate install needed for them.

For very large indexes (whole genome with introns, atlas-scale references), you may need to build kallisto from source for the v0.50+ features:

git clone https://github.com/pachterlab/kallisto && cd kallisto && mkdir build && cd build
cmake .. -DUSE_HDF5=ON && make && make install

Pipeline 1 — Build the Reference Index

You only do this once per genome. Reuse for all your samples.

Standard transcriptome (most users)

# Inputs: genome FASTA + GTF (matching versions)
# Example with GENCODE human v44:
kb ref \
    -i  index.idx \
    -g  t2g.txt \
    -f1 cdna.fasta \
    GRCh38.primary_assembly.genome.fa.gz \
    gencode.v44.annotation.gtf.gz

Output files:

• index.idx — the kallisto pseudoalignment index (~3-5 GB for human)
• t2g.txt — transcript-to-gene mapping
• cdna.fasta — extracted transcript sequences (intermediate, keep if rebuilding)

This takes ~30-60 minutes for human/mouse on a typical workstation. Index files are reusable across all your samples.

Pre-built indexes (skip the build)

For common references, kb-python can download pre-built indexes:

kb ref -d human  -i index.idx -g t2g.txt
kb ref -d mouse  -i index.idx -g t2g.txt
kb ref -d zebrafish -i index.idx -g t2g.txt

Faster than building, identical result. Recommended unless you have a custom genome.

nac workflow — nascent + mature RNA (for snRNA-seq and RNA velocity)

Single-nucleus RNA-seq captures predominantly intronic reads (from pre-mRNA in the nucleus). Standard transcriptome indexes only handle exonic/mature mRNA, so you lose those reads. The nac workflow indexes both.

kb ref \
    --workflow nac \
    -i  index.idx \
    -g  t2g.txt \
    -c1 cdna.txt          \
    -c2 nascent.txt       \
    -f1 cdna.fasta        \
    -f2 nascent.fasta     \
    GRCh38.primary_assembly.genome.fa.gz \
    gencode.v44.annotation.gtf.gz

Additional outputs vs standard:

• cdna.txt — transcript IDs for mature (spliced) sequences
• nascent.txt — transcript IDs for nascent (intron-containing) sequences
• nascent.fasta — pre-mRNA sequences for the index

Required when:

• Running snRNA-seq (most nuclei reads land in introns)
• Doing RNA velocity (need separate spliced + unspliced counts)

For RNA velocity specifically, also consider the older lamanno workflow which produces the spliced/unspliced layout scVelo expects:

kb ref --workflow lamanno -i index.idx -g t2g.txt \
    -f1 cdna.fasta -f2 intron.fasta \
    -c1 cdna_t2c.txt -c2 intron_t2c.txt \
    GRCh38.primary_assembly.genome.fa.gz gencode.v44.annotation.gtf.gz

Convenience: bash scripts/kb_ref.sh --genome GRCh38.fa --gtf annotation.gtf --out-dir reference/. See references/advanced.md for nac vs lamanno trade-offs and how to choose.

Source: Index generation docs.

Pipeline 2 — Quantify a Sample

# Standard 10X v3 (cellranger output style)
kb count \
    -i index.idx \
    -g t2g.txt \
    -x 10xv3 \
    -o output_dir \
    --h5ad \
    sample_R1.fastq.gz sample_R2.fastq.gz

Key arguments

Supported chemistries (-x values)

Full list: kb --list after installation. For unsupported chemistries, define one in seqspec and pass via -x <seqspec.yaml>.

Input FASTQ layout

For paired-end droplet chemistries, file order is R1 first, R2 second, repeated per lane:

kb count -i index.idx -g t2g.txt -x 10xv3 -o out --h5ad \
    L001_R1.fastq.gz L001_R2.fastq.gz \
    L002_R1.fastq.gz L002_R2.fastq.gz \
    L003_R1.fastq.gz L003_R2.fastq.gz

For 10X v1 (legacy), also include the I1 file: I1 R1 R2 I1 R1 R2 ....

For SMART-Seq2 / full-length protocols where both reads are biological, use --parity paired:

kb count -i index.idx -g t2g.txt -x SMARTSEQ2 -o out --parity paired \
    SAMPLE1_1.fastq.gz SAMPLE1_2.fastq.gz \
    SAMPLE2_1.fastq.gz SAMPLE2_2.fastq.gz

Output structure

output_dir/
├── counts_unfiltered/                 # cell × gene matrix BEFORE cell calling
│   ├── cells_x_genes.mtx               # sparse matrix (Matrix Market format)
│   ├── cells_x_genes.barcodes.txt     # cell barcodes (one per row)
│   ├── cells_x_genes.genes.txt        # gene IDs (Ensembl)
│   ├── cells_x_genes.genes.names.txt  # gene symbols
│   └── adata.h5ad                     # AnnData (if --h5ad was passed)
├── inspect.json                        # QC metrics from bustools
├── run_info.json                       # kallisto run info (alignment rate, etc.)
├── output.bus                          # raw BUS file
├── output.unfiltered.bus               # barcode-corrected
└── matrix.ec                           # equivalence classes

The counts_unfiltered/ directory is your handoff point. Empty droplets are NOT filtered out at this stage — call cells with DropletUtils::emptyDrops or CellBender downstream.

Convenience: bash scripts/kb_count.sh --index index.idx --t2g t2g.txt --tech 10xv3 --out-dir sample1_out fastq/R1.fastq.gz fastq/R2.fastq.gz.

Source: pseudoalignment + output files docs.

Pipeline 3 — Handoff to Downstream (scanpy / Seurat / SCE)

scanpy / Python

import scanpy as sc
import anndata

# Easiest: load the .h5ad written by kb count --h5ad
adata = anndata.read_h5ad('output_dir/counts_unfiltered/adata.h5ad')

# Or load the .mtx manually if --h5ad wasn't passed
adata = sc.read_mtx('output_dir/counts_unfiltered/cells_x_genes.mtx')
barcodes = open('output_dir/counts_unfiltered/cells_x_genes.barcodes.txt').read().splitlines()
genes    = open('output_dir/counts_unfiltered/cells_x_genes.genes.txt').read().splitlines()
gene_names = open('output_dir/counts_unfiltered/cells_x_genes.genes.names.txt').read().splitlines()
adata.obs_names = barcodes
adata.var_names = gene_names
adata.var['gene_id'] = genes

# Standard preprocessing (see scanpy protocol for the full pipeline)
sc.pp.filter_cells(adata, min_genes=200)
sc.pp.filter_genes(adata, min_cells=3)
sc.pp.normalize_total(adata, target_sum=1e4)
sc.pp.log1p(adata)

Seurat / R

library(Seurat)
library(Matrix)

mtx <- readMM("output_dir/counts_unfiltered/cells_x_genes.mtx")
barcodes <- readLines("output_dir/counts_unfiltered/cells_x_genes.barcodes.txt")
genes    <- readLines("output_dir/counts_unfiltered/cells_x_genes.genes.names.txt")

# kb writes cells × genes; Seurat expects genes × cells — transpose
mtx <- t(mtx)
rownames(mtx) <- genes
colnames(mtx) <- barcodes
mtx <- as(mtx, "CsparseMatrix")

obj <- CreateSeuratObject(counts = mtx, min.cells = 3, min.features = 200)
# Now standard Seurat from here — see seurat protocol

SingleCellExperiment / R

library(SingleCellExperiment)
library(Matrix)

mtx <- t(readMM("output_dir/counts_unfiltered/cells_x_genes.mtx"))
rownames(mtx) <- readLines("output_dir/counts_unfiltered/cells_x_genes.genes.names.txt")
colnames(mtx) <- readLines("output_dir/counts_unfiltered/cells_x_genes.barcodes.txt")

sce <- SingleCellExperiment(assays = list(counts = mtx))
# Now standard SCE pipeline — see singlecellexperiment protocol

Source: kb → scanpy notebook.

nac Workflow — Loading Nascent + Mature Counts

When you ran kb count --workflow nac, the output has BOTH a mature (cdna) and nascent matrix:

import scanpy as sc

# Each matrix is loaded separately
mature  = sc.read_mtx('output_dir/counts_unfiltered/cells_x_genes.mature.mtx')
nascent = sc.read_mtx('output_dir/counts_unfiltered/cells_x_genes.nascent.mtx')

# For RNA velocity, stash both as layers of one AnnData
import anndata
adata = anndata.AnnData(
    X      = mature.X,
    layers = {'spliced': mature.X, 'unspliced': nascent.X},
    obs    = mature.obs,
    var    = mature.var,
)
# Then pass to scVelo (see scvelo protocol)

For snRNA-seq (no velocity, just better gene counts), sum spliced + unspliced before downstream — nuclei reads come from both.

Multi-Sample Batch Processing

# Create a batch file: SAMPLE_ID<TAB>R1_PATH<TAB>R2_PATH per row
cat > samples.txt <<EOF
sample1	/data/s1_R1.fastq.gz	/data/s1_R2.fastq.gz
sample2	/data/s2_R1.fastq.gz	/data/s2_R2.fastq.gz
sample3	/data/s3_R1.fastq.gz	/data/s3_R2.fastq.gz
EOF

# Run all samples in one kb count call — preserves sample identity
kb count -i index.idx -g t2g.txt -x 10xv3 -o multi_out --h5ad \
    --batch-barcodes \
    --tcc \
    batch.txt

Each cell's barcode gets prefixed with its sample ID, so you can split later.

Key Parameters to Adjust

kb ref

• --workflow — standard (default), nac (nascent+mature for snRNA-seq + velocity), lamanno (legacy velocity layout)
• Memory: index build for human/mouse needs ~16 GB; use a smaller k-mer (default 31) only for tiny references

kb count

• -x <tech> — must match your library chemistry exactly; use kb --list to check
• --workflow — must match the kb ref workflow
• --h5ad — recommended; writes a ready-to-load AnnData
• --filter bustools — applies bustools knee-plot cell calling (optional)
• -t — threads; one per core, kb scales near-linearly to ~8 threads
• -m 8G — RAM cap; bump if you see OOM errors

Memory considerations

• Index files for human/mouse: ~3-5 GB. Keep them on local SSD for kb count.
• kb count peak RAM is roughly 1.5× the index size.
• For atlas-scale runs (1000+ samples), build the index once and reuse across all kb count jobs.

Best Practices

• Build the index once, reuse forever. It's deterministic — same FASTA + GTF + workflow gives the same index. Stash it on shared storage.
• Use --h5ad. Saves writing the awkward MTX → AnnData glue code. Works for scanpy and converts cleanly to Seurat / SCE.
• Don't filter cells in kb count. Skip --filter bustools for most analyses — use scanpy min_genes=200 + min_cells=3 downstream, or CellBender for ambient-aware filtering.
• For snRNA-seq, ALWAYS use --workflow nac. Without nascent indexing, you throw away 30-70% of reads (the introns).
• Check run_info.json after every run. Pseudoalignment rate < 70% means your reference doesn't match the data well — check species, gene model version.
• Multi-sample: use --batch-barcodes + a samples file. One kb count call across many samples is much faster than parallel jobs because the index loads once.
• Reproducibility: pin kb-python + kallisto + bustools versions with kb --version. Index files can subtly change across versions.

Common Pitfalls

• Mismatched kb ref and kb count workflows. kb ref --workflow nac + kb count without --workflow nac silently produces wrong results.
• Forgetting -x chemistry. Default is 10x v3; using it on v2 / Drop-seq data fails silently (most barcodes won't correct).
• Confusing R1 (barcode + UMI) with R2 (cDNA). Most chemistries are R1=barcode, R2=cDNA. The order in the command matters.
• Using filtered_feature_bc_matrix output style as input to kb. kb takes raw FASTQs; it's the upstream step. cellranger's matrix outputs aren't compatible.
• Index built on wrong gene model. If counts look 30% lower than expected from cellranger, the GTF probably annotates fewer transcripts. Use the same GTF cellranger uses (typically GENCODE or Ensembl primary assembly) for fair comparison.

End-to-End Template

assets/kallisto_template.sh — single parameterized script. Edit FASTQ paths, technology, and reference paths, then run end-to-end ref → count → handoff message.

Convenience Scripts

• scripts/kb_ref.sh — index builder (standard or nac workflow)
• scripts/kb_count.sh — per-sample quantifier

For downstream analysis, switch to one of:

• scanpy — Python downstream
• seurat — R / Seurat downstream
• singlecellexperiment — R / Bioconductor downstream
• cellbender — ambient-aware cell calling on the unfiltered matrix
• scvelo — RNA velocity if you used --workflow nac or lamanno

References

• kallisto docs — Pachter Lab
• Introduction
• Installation
• Index generation
• Single-cell quantification
• scanpy notebook tutorial
• Bray et al. (2016), Near-optimal probabilistic RNA-seq quantification, Nature Biotechnology (original kallisto paper)
• Melsted et al. (2021), Modular, efficient and constant-memory single-cell RNA-seq preprocessing, Nature Biotechnology (kb-python paper)
• Sullivan et al. (2024), kallisto, bustools, and kb-python for quantifying bulk, single-cell, and single-nucleus RNA-seq, Nature Protocols

Source: swaruplab/operon — distributed by TomeVault.