Quick-start tutorial

This section contains quick start tutorials for different CellBender modules.

remove-background

In this tutorial, we will run remove-background on a small dataset derived from the 10x Genomics heart10k snRNA-seq dataset (v3 Chemistry, CellRanger 3.0.2).

As a first step, we download the full dataset and generate a smaller trimmed copy by selecting 500 barcodes with high UMI count (likely non-empty) and an additional 50,000 barcodes with small UMI count (likely empty). We also trim to keep only the top 100 most highly-expressed genes. Note that the trimming step is performed in order to allow us go through this tutorial in a minute on a typical CPU. Processing the full untrimmed dataset requires a CUDA-enabled GPU (e.g. NVIDIA Testla T4) and takes about 30 minutes to finish.

(Please note that trimming is NOT part of the recommended workflow, and is only for the purposes of a quick demo run. When you run remove-background on real data, do not do any preprocessing or feature selection or barcode selection first.)

We have created a python script to download and trim the dataset. Navigate to examples/remove_background/ under your CellBender installation root directory and run the following command in the console:

$ python generate_tiny_10x_dataset.py

After successful completion of the script, you should have a new file named tiny_raw_feature_bc_matrix.h5ad.

Run remove-background

We proceed to run remove-background on the trimmed dataset using the following command:

(cellbender) $ cellbender remove-background \
     --input tiny_raw_feature_bc_matrix.h5ad \
     --output tiny_output.h5 \
     --expected-cells 500 \
     --total-droplets-included 2000

Again, here we leave out the --cuda flag solely for the purposes of being able to run this command on a CPU. But a GPU is highly recommended for real datasets.

The computation will finish within a minute or two (after 150 epochs). The tool outputs the following files:

tiny_output.h5: An HDF5 file containing a detailed output of the inference procedure, including the normalized abundance of ambient transcripts, contamination fraction of each droplet, a low-dimensional embedding of the background-corrected gene expression, and the background-corrected counts matrix (in CSC sparse format). Please refer to the full documentation for a detailed description of these and other fields.
tiny_output_filtered.h5: Same as above, though, only including droplets with a posterior cell probability exceeding 0.5.
tiny_output_cell_barcodes.csv: The list of barcodes with a posterior cell probability exceeding 0.5.
tiny_output.pdf: A PDF summary of the results showing (1) the evolution of the loss function during training, (2) a ranked-ordered total UMI plot along with posterior cell probabilities, and (3) a two-dimensional PCA scatter plot of the latent embedding of the expressions in cell-containing droplets. Notice the rapid drop in the cell probability after UMI rank ~ 500.
tiny_output_umi_counts.pdf: A PDF showing the UMI counts per droplet as a histogram, annotated with what CellBender thinks is empty versus cells. This describes CellBender’s prior. This is mainly a diagnostic if something seems to be going wrong with CellBender’s automatic prior determination.
tiny_output.log: Log file.
tiny_output_metrics.csv: Output metrics, most of which have very descriptive names. This file is not used by most users, but the idea is that it can be incorporated into automated pipelines which could re-run CellBender automatically (with different parameters) if something goes wrong.
tiny_output_report.html: An HTML report which points out a few things about the run and highlights differences between the output and the input. Issues warnings if there are any aspects of the run that look anomalous, and makes suggestions.

Finally, try running the tool with --expected-cells 100 and --expected-cells 1000. You should find that the output remains virtually the same.

Use output count matrix in downstream analyses

The count matrix that remove-background generates can be easily loaded and used for downstream analyses in scanpy and Seurat.

You can load the filtered count matrix (containing only cells) into scanpy:

# import scanpy
import scanpy as sc

# load the data
adata = sc.read_10x_h5('tiny_output_filtered.h5')

This will yield adata like this

AnnData object with n_obs × n_vars = 531 × 100
    var: 'gene_ids', 'feature_types', 'genome'

The CellBender output counts are in adata.X.

However, this does not include any of the CellBender metadata when loading the file. To include the metadata, but still load an AnnData object that scanpy can operate on, try some of the functions from cellbender.remove_background.downstream (see here)

# import function
from cellbender.remove_background.downstream import anndata_from_h5

# load the data
adata = anndata_from_h5('tiny_output.h5')

This yields an adata with all the cell barcodes which were analyzed by CellBender (all the --total-droplets-included), along with all the metadata and latent variables inferred by CellBender:

AnnData object with n_obs × n_vars = 1000 × 100
    obs: 'background_fraction', 'cell_probability', 'cell_size', 'droplet_efficiency'
    var: 'ambient_expression', 'features_analyzed_inds', 'feature_type', 'genome', 'gene_id'
    uns: 'cell_size_lognormal_std', 'empty_droplet_size_lognormal_loc', 'empty_droplet_size_lognormal_scale', 'posterior_regularization_lambda', 'swapping_fraction_dist_params', 'target_false_positive_rate', 'fraction_data_used_for_testing', 'test_elbo', 'test_epoch', 'train_elbo', 'train_epoch'
    obsm: 'gene_expression_encoding'

(If you want to load both the raw data and the CellBender data into one AnnData object, which is very useful, try the load_anndata_from_input_and_output() function in cellbender.remove_background.downstream, see here)

You can access the latent gene expression embedding learned by CellBender in adata.obsm['gene_expression_encoding'], the inferred ambient RNA profile is in adata.var['ambient_expression'], and the inferred cell probabilties are in adata.obs['cell_probability'].

You can limit this adata to CellBender cell calls very easily:

adata[adata.obs['cell_probability'] > 0.5]

View of AnnData object with n_obs × n_vars = 531 × 100
    obs: 'background_fraction', 'cell_probability', 'cell_size', 'droplet_efficiency'
    var: 'ambient_expression', 'features_analyzed_inds', 'feature_type', 'genome', 'gene_id'
    uns: 'cell_size_lognormal_std', 'empty_droplet_size_lognormal_loc', 'empty_droplet_size_lognormal_scale', 'posterior_regularization_lambda', 'swapping_fraction_dist_params', 'target_false_positive_rate', 'fraction_data_used_for_testing', 'test_elbo', 'test_epoch', 'train_elbo', 'train_epoch'
    obsm: 'gene_expression_encoding'

How to use the latent gene expression downstream

After loading data using the anndata_from_h5() function as shown above, we can compute nearest neighbors in scanpy, using the CellBender latent representation of cells, and make a UMAP and do clustering:

# compute a UMAP and do clustering using the cellbender latent gene expression embedding
sc.pp.neighbors(adata, use_rep='gene_expression_encoding', metric='euclidean', method='umap')
sc.pp.umap(adata)
sc.pp.leiden(adata)

Seurat

Seurat 4.0.2 uses a dataloader Read10X_h5() which is not currently compatible with the CellBender output file format. Hopefully Seurat will update its dataloader to ignore extra information in the future, but in the interim, we can use a super handy utility from PyTables to strip the extra CellBender information out of the output file so that Seurat can load it.

From a python environment in which PyTables is installed, do the following at the command line:

$ ptrepack --complevel 5 tiny_output_filtered.h5:/matrix tiny_output_filtered_seurat.h5:/matrix

(The flag --complevel 5 ensures that the file size does not increase.)

The file tiny_output_filtered_seurat.h5 is now formatted exactly like a CellRanger v3 h5 file, so Seurat can load it:

# load data from the filtered h5 file
data.file <- 'tiny_output_filtered_seurat.h5'
data.data <- Read10X_h5(filename = data.file, use.names = TRUE)

# create Seurat object
obj <- CreateSeuratObject(counts = data.data)
obj

An object of class Seurat
100 features across 531 samples within 1 assay
Active assay: RNA (100 features, 0 variable features)

Of course, this will not load any metadata from CellBender, so if that is desired, it would have to be accessed and added to the object another way.

Another option for loading data into Seurat would be third party packages like scCustomize from Samuel Marsh.