kb-python

kb-python is a python package for processing single-cell RNA-sequencing. It wraps the kallisto | bustools single-cell RNA-seq command line tools in order to unify multiple processing workflows.

kb-python was first developed by Kyung Hoi (Joseph) Min and A. Sina Booeshaghi while in Lior Pachter's lab at Caltech. If you use kb-python in a publication please cite*:

Melsted, P., Booeshaghi, A.S., et al. 
Modular, efficient and constant-memory single-cell RNA-seq preprocessing. 
Nat Biotechnol  39, 813–818 (2021). 
https://doi.org/10.1038/s41587-021-00870-2

Installation

The latest release can be installed with

pip install kb-python

The development version can be installed with

pip install git+https://github.com/pachterlab/kb_python

There are no prerequisite packages to install. The kallisto and bustools binaries are included with the package.

Usage

kb consists of six subcommands

$ kb
usage: kb [-h] [--list] <CMD> ...
positional arguments:
  <CMD>
    info      Display package and citation information
    compile   Compile `kallisto` and `bustools` binaries from source
    ref       Build a kallisto index and transcript-to-gene mapping
    count     Generate count matrices from a set of single-cell FASTQ files
    extract   Extract reads that were pseudoaligned to specific genes/transcripts (or extract all reads that were / were not pseudoaligned)
    sweep     Remove ambient counts from the count matrix with cellsweep

kb ref: generate a pseudoalignment index

The kb ref command takes in a species annotation file (GTF) and associated genome (FASTA) and builds a species-specific index for pseudoalignment of reads. This must be run before kb count. Internally, kb ref extracts the coding regions from the GTF and builds a transcriptome FASTA that is then indexed with kallisto index.

kb ref -i index.idx -g t2g.txt -f1 transcriptome.fa <GENOME> <GENOME_ANNOTATION>

• <GENOME> refers to a genome file (FASTA).
  • For example, the zebrafish genome is hosted by ensembl and can be downloaded here
• <GENOME_ANNOTATION> refers to a genome annotation file (GTF)
  • For example, the zebrafish genome annotation file is hosted by ensembl and can be downloaded here
• Note: The latest genome annotation and genome file for every species on ensembl can be found with the gget command-line tool.

Prebuilt indices are available at https://github.com/pachterlab/kallisto-transcriptome-indices

Examples

# Index the transcriptome from genome FASTA (genome.fa.gz) and GTF (annotation.gtf.gz)
$ kb ref -i index.idx -g t2g.txt -f1 transcriptome.fa genome.fa.gz annotation.gtf.gz
# An example for downloading a prebuilt reference for mouse
$ kb ref -d mouse -i index.idx -g t2g.txt

---

kb count: pseudoalign and count reads

The kb count command takes in the pseudoalignment index (built with kb ref) and sequencing reads generated by a sequencing machine to generate a count matrix. Internally, kb count runs numerous kallisto and bustools commands comprising a single-cell workflow for the specified technology that generated the sequencing reads.

kb  count -i index.idx -g t2g.txt -o out/ -x <TECHNOLOGY> <FASTQ FILE[s]>

• <TECHNOLOGY> refers to the assay that generated the sequencing reads.
  • For a list of supported assays run kb --list
• <FASTQ FILE[s]> refers to the a list of FASTQ files generated
  • Different assays will have a different number of FASTQ files
  • Different assays will place the different features in different FASTQ files
    • For example, sequencing a 10xv3 library on a NextSeq Illumina sequencer usually results in two FASTQ files.
    • The R1.fastq.gz file (colloquially called "read 1") contains a 16 basepair cell barcode and a 12 basepair unique molecular identifier (UMI).
    • The R2.fastq.gz file (colloquially called "read 2") contains the cDNA associated with the cell barcode-UMI pair in read 1.

Examples

# Quantify 10xv3 reads read1.fastq.gz and read2.fastq.gz
$ kb count -i index.idx -g t2g.txt -o out/ -x 10xv3 read1.fastq.gz read2.fastq.gz

---

kb info: display package and citation information

The kb info command prints out package information including the version of kb-python, kallisto, and bustools along with their installation location.

$ kb info
kb_python 0.31.0 ...
kallisto: 0.52.0 ...
bustools: 0.45.1 ...
...

---

kb compile: compile kallisto and bustools binaries from source

The kb compile command grabs the latest kallisto and bustools source and compiles the binaries. Note: this is not required to run kb-python.

---

kb sweep: remove ambient RNA with cellsweep

The kb sweep command removes ambient and bulk RNA contamination from the unfiltered count matrix produced by kb count, using cellsweep. cellsweep is not installed with kb-python; install the required version with pip install cellsweep[analysis]==1.0.2. See How kb sweep works for details.

kb sweep [-c CELLTYPES] [--h5ad] [--celltype-column COLUMN] [--leiden-resolution RESOLUTION] -o OUT <KB_COUNT_DIR>

• <KB_COUNT_DIR> refers to the output directory of kb count (or an .h5ad file with the unfiltered count matrix)
• --h5ad reads counts_unfiltered/adata.h5ad instead of the .mtx and associated files
• -c is a tab-delimited file containing barcodecelltype, with an optional third column giving the empty status (true/false). Barcodes in the count matrix but not in the file are treated as empty droplets
• --celltype-column is the column in the h5ad obs holding celltypes (default: celltype)
• --leiden-resolution assigns celltypes by Leiden clustering at the given resolution
• --keep-empties keeps the empty droplets in the output (by default only the real cells are written)
• --round-X rounds the denoised counts to integers
• kb sweep --help lists the empty-droplet options

Examples

# Denoise using a barcode-to-celltype mapping
$ kb sweep -c celltypes.txt -o out/adata_denoised.h5ad out/
# Denoise the h5ad from kb count --h5ad, using its "celltype" column
$ kb sweep --h5ad out/
# Denoise, assigning celltypes by Leiden clustering
$ kb sweep --leiden-resolution 1.0 out/

Use cases

kb-python facilitates fast and uniform pre-processing of single-cell sequencing data to answer relevant research questions.

$ pip install kb-python gget ffq

# Goal: quantify publicly available scRNAseq data
$ kb ref -i index.idx -g t2g.txt -f1 transcriptome.fa $(gget ref --ftp -w dna,gtf homo_sapiens)
$ kb count -i index.idx -g t2g.txt -x 10xv3 -o out $(ffq --ftp SRR10668798 | jq -r '.[] | .url' | tr '\n' ' ')
# -> count matrix in out/ folder

# Goal: quantify 10xv2 feature barcode data, feature_barcodes.txt is a tab-delimited file
# containing barcode_sequence<tab>barcode_name
$ kb ref -i index.idx -g f2g.txt -f1 features.fa --workflow kite feature_barcodes.txt
$ kb count -i index.idx -g f2b.txt -x 10xv2 -o out/ --workflow kite --h5ad R1.fastq.gz R2.fastq.gz
# -> count matrix in out/ folder

Submitted by @sbooeshaghi.

Do you have a cool use case for kb-python? Submit a PR (including the goal, code snippet, and your username) so that we can feature it here.

Documentation and Tutorials

Documentation and tutorials for kb-python usage are available at https://kallisto.readthedocs.io/en/latest/.

Contributing

Thank you for wanting to improve kb-python! If you have believe you've found a bug, please submit an issue.

If you have a new feature you'd like to add to kb-python please create a pull request. Pull requests should contain a message detailing the exact changes made, the reasons for the change, and tests that check for the correctness of those changes.

Cite

If you use kb-python in a publication, please cite the following papers:

kb-python & kallisto and/or bustools

@article{sullivan2023kallisto,
  title={kallisto, bustools, and kb-python for quantifying bulk, single-cell, and single-nucleus RNA-seq},
  author={Sullivan, Delaney K and Min, Kyung Hoi and Hj{\"o}rleifsson, Kristj{\'a}n Eldj{\'a}rn and Luebbert, Laura and Holley, Guillaume and Moses, Lambda and Gustafsson, Johan and Bray, Nicolas L and Pimentel, Harold and Booeshaghi, A Sina and others},
  journal={bioRxiv},
  pages={2023--11},
  year={2023},
  publisher={Cold Spring Harbor Laboratory}
}

bustools

@article{melsted2021modular,
  title={\href{https://doi.org/10.1038/s41587-021-00870-2}{Modular, efficient and constant-memory single-cell RNA-seq preprocessing}},
  author={Melsted, P{\'a}ll and Booeshaghi, A. Sina and Liu, Lauren and Gao, Fan and Lu, Lambda and Min, Kyung Hoi Joseph and da Veiga Beltrame, Eduardo and Hj{\"o}rleifsson, Kristj{\'a}n Eldj{\'a}rn and Gehring, Jase and Pachter, Lior},
  author+an={1=first;2=first,highlight},
  journal={Nature biotechnology},
  year={2021},
  month={4},
  day={1},
  doi={https://doi.org/10.1038/s41587-021-00870-2}
}

kallisto

@article{bray2016near,
  title={Near-optimal probabilistic RNA-seq quantification},
  author={Bray, Nicolas L and Pimentel, Harold and Melsted, P{\'a}ll and Pachter, Lior},
  journal={Nature biotechnology},
  volume={34},
  number={5},
  pages={525--527},
  year={2016},
  publisher={Nature Publishing Group}
}

kITE

@article{booeshaghi2024quantifying,
  title={Quantifying orthogonal barcodes for sequence census assays},
  author={Booeshaghi, A Sina and Min, Kyung Hoi and Gehring, Jase and Pachter, Lior},
  journal={Bioinformatics Advances},
  volume={4},
  number={1},
  pages={vbad181},
  year={2024},
  publisher={Oxford University Press}
}

BUS format

@article{melsted2019barcode,
  title={The barcode, UMI, set format and BUStools},
  author={Melsted, P{\'a}ll and Ntranos, Vasilis and Pachter, Lior},
  journal={Bioinformatics},
  volume={35},
  number={21},
  pages={4472--4473},
  year={2019},
  publisher={Oxford University Press}
}

kb-python was inspired by Sten Linnarsson’s loompy fromfq command (http://linnarssonlab.org/loompy/kallisto/index.html)

How kb sweep works

kb sweep is a wrapper around cellsweep.denoise_count_matrix. It runs in three steps.

1. Load the count matrix. If the input is a kb count output directory, kb sweep reads the unfiltered matrix from counts_unfiltered/: adata.h5ad if --h5ad is given, otherwise cells_x_genes.mtx with cells_x_genes.barcodes.txt and cells_x_genes.genes.names.txt. An .h5ad file can also be given directly, and is treated as h5ad input. By default the output is written to counts_unfiltered/adata_denoised.h5ad (or <input>_denoised.h5ad for an .h5ad file); use --overwrite to replace an existing output.

2. Assign celltypes. cellsweep needs a celltype for each cell, which it reads from a column of adata.obs named by its celltype_key argument. kb sweep takes the celltypes from the first of these sources that is available:

1. The -c file, if provided. Each line holds a barcode, its celltype and, optionally, its empty status (true/false, also t/f, yes/no, 1/0), separated by tabs. Either every line or no line has the third column, and the celltype may be blank only for a barcode marked empty. The celltypes are stored in adata.obs["celltype"] (replacing any existing column of that name). Barcodes in the count matrix but not in the file are treated as empty droplets, and barcodes in the file but not in the count matrix are ignored; a warning reports how many of each, and it is an error if no barcodes match. The resulting empty status is stored in adata.obs["is_empty"], which takes the place of empty droplet detection (--empty-droplet-method, --umi-cutoff and --expected-cells are then unused). If the file has no third column and lists every barcode, it says nothing about which are empty, so empty droplets are detected as usual.
2. The --celltype-column column (default celltype) of the input AnnData's obs. This is only checked for h5ad input, since the .mtx files carry no celltype annotations. The column is passed to cellsweep as-is; it is not copied.
3. Leiden clustering, if --leiden-resolution is provided. A copy of the matrix is preprocessed and clustered with Scanpy (2,000 highly variable genes, 50 principal components, 15 neighbors), and the Leiden clusters are stored in adata.obs["celltype"]. The raw counts that are denoised are left untouched.

If none of these is available, kb sweep exits with an error asking for -c, a celltype column, or --leiden-resolution.

3. Run cellsweep. kb sweep first checks that cellsweep version 1.0.2 is installed. It then calls cellsweep.denoise_count_matrix. This identifies empty droplets (with --empty-droplet-method, --umi-cutoff and --expected-cells, unless adata.obs["is_empty"] is already set) and computes a mean expression profile for each celltype. It then fits an Expectation-Maximization model in which each observed count is a mixture of ambient RNA, bulk RNA, and the cell's true celltype signal. The model estimates the ambient fraction of each cell, the bulk contamination fraction, the celltype profiles and the ambient profile, and it stops when the log-likelihood, the celltype profiles and the contamination fractions have all stabilized or after a maximum number of iterations. The denoised matrix is stored in adata.X, with the raw counts in adata.layers["raw"], the per-cell ambient fraction (alpha_hat) and total contamination fraction (contamination_fraction) in adata.obs, the fitted ambient and bulk profiles in adata.var, and the celltype profiles in adata.uns. By default the empty droplets are removed from the output after fitting; --keep-empties retains them, and --round-X rounds the denoised counts to integers.