GPU-Accelerated Single-Cell Genomics Analysis with RAPIDS

This repository contains example notebooks demonstrating how to use RAPIDS for GPU-accelerated analysis of single-cell sequencing data.

RAPIDS is a suite of open-source Python libraries that can speed up data science workflows using GPU acceleration.Starting from a single-cell count matrix, RAPIDS libraries can be used to perform data processing, dimensionality reduction, clustering, visualization, and comparison of cell clusters.

Several of our examples are inspired by the Scanpy tutorials and based upon the AnnData format. Currently, we provide examples for scRNA-seq and scATAC-seq, and we have scaled up to 1 million cells. We also show how to create GPU-powered interactive, in-browser visualizations to explore single-cell datasets.

Dataset sizes for single-cell genomics studies are increasing, presently reaching millions of cells. With RAPIDS, it becomes easy to analyze large datasets interactively and in real time, enabling faster scientific discoveries.

Installation

Docker container

A container with all dependencies, notebooks and source code is available at https://hub.docker.com/r/claraparabricks/single-cell-examples_rapids_cuda11.0.

Please execute the following commands to start the notebook and follow the URL in the log to open Jupyter web application.

docker pull claraparabricks/single-cell-examples_rapids_cuda11.0

docker run --gpus all --rm -v /mnt/data:/data claraparabricks/single-cell-examples_rapids_cuda11.0

conda

All dependencies for these examples can be installed with conda.

conda env create --name rapidgenomics -f conda/rapidgenomics_cuda11.5.yml
conda activate rapidgenomics
python -m ipykernel install --user --display-name "Python (rapidgenomics)"

After installing the necessary dependencies, you can just run jupyter lab. There are are a few different conda environment files which correspond to different notebooks. In addition to the one listed above, there is one for the CPU notebooks, one for the real-time visualization notebook, and one for the AtacSeq notebook.

Configuration

Unified Virtual Memory (UVM) can be used to oversubscribe your GPU memory so that chunks of data will be automatically offloaded to main memory when necessary. This is a great way to explore data without having to worry about out of memory errors, but it does degrade performance in proportion to the amount of oversubscription. UVM is enabled by default in these examples and can be enabled/disabled in any RAPIDS workflow with the following:

import cupy as cp
import rmm
rmm.reinitialize(managed_memory=True)
cp.cuda.set_allocator(rmm.rmm_cupy_allocator)

RAPIDS provides a GPU Dashboard, which contains useful tools to monitor GPU hardware right in Jupyter.

Citation

If you use this code, please cite our preprint:

Nolet C., Lal A., et al., (2022). Accelerating single-cell genomic analysis with GPUs. bioRxiv.

Example 1: Single-cell RNA-seq of 70,000 Human Lung Cells

We use RAPIDS to accelerate the analysis of a ~70,000-cell single-cell RNA sequencing dataset from human lung cells. This example includes preprocessing, dimension reduction, clustering, visualization and gene ranking.

Example Dataset

The dataset is from Travaglini et al. 2020. If you wish to run the example notebook using the same data, use the following command to download the count matrix for this dataset and store it in the data folder:

wget -P <path to this repository>/data https://rapids-single-cell-examples.s3.us-east-2.amazonaws.com/krasnow_hlca_10x.sparse.h5ad

Example Code

Follow this Jupyter notebook for RAPIDS analysis of this dataset. In order for the notebook to run, the file rapids_scanpy_funcs.py needs to be in the same folder as the notebook.

We provide a second notebook with the CPU version of this analysis here.

Acceleration

We report the runtime of these notebooks on various GCP instances below. All runtimes are given in seconds. Acceleration is given in parentheses. Benchmarking was performed on Dec 16, 2020.

Step	CPU n1-standard-16 16 vCPUs	GPU n1-standard-16 T4 16 GB GPU (Acceleration)	GPU n1-highmem-8 Tesla V100 16 GB GPU (Acceleration)	GPU a2-highgpu-1g Tesla A100 40GB GPU (Acceleration)
Preprocessing	70	92 (0.76x)	53 (1.3x)	59 (1.2x)
PCA	10.6	5.0 (2.1x)	3.2 (3.3x)	2.7 (3.9x)
t-SNE	220	2.8 (79x)	1.4 (157x)	2.2 (100x)
k-means (single iteration)	14.3	0.31 (46x)	0.12 (119x)	0.08 (179x)
KNN	17.8	7.6 (2.3x)	6.8 (2.6x)	5.7 (3.1x)
UMAP	97	1.1 (88x)	0.55 (176x)	0.53 (x)
Louvain clustering	13.9	0.24 (58x)	0.15 (93x)	0.11 (126x)
Leiden clustering	12.8	0.17 (75x)	0.09 (142x)	0.08 (160x)
Differential Gene Expression	153	36 (4.3x)	7.5 (20.4x)	6.3 (24x)
Re-analysis of subgroup	29	7.6 (3.8x)	3.6 (8x)	3.5 (8x)
End-to-end notebook run	654	166	93	96
Price ($/hr)	0.760	1.110	2.953	3.673
Total cost ($)	0.138	0.051	0.076	0.098

Example 2: Single-cell RNA-seq of 1.3 Million Mouse Brain Cells

We demonstrate the use of RAPIDS to accelerate the analysis of single-cell RNA-seq data from 1.3 million cells. This example includes preprocessing, dimension reduction, clustering and visualization.

Compared to the previous example, here we make several adjustments to handle the larger dataset. We perform most of the preprocessing operations (e.g. filtering, normalization) while reading the dataset in batches. Further, we perform a batched PCA by training on a fraction of cells and transforming the data in batches.

This example relies heavily on UVM. While it should work on any GPU built on the Pascal architecture or newer, you will want to make sure there is enough main memory available. Oversubscribing a GPU by more than a factor of 2x can cause thrashing in UVM, which can ultimately lead to the notebook freezing.

Example Dataset

The dataset was made publicly available by 10X Genomics. Use the following command to download the count matrix for this dataset and store it in the data folder:

wget -P <path to this repository>/data https://rapids-single-cell-examples.s3.us-east-2.amazonaws.com/1M_brain_cells_10X.sparse.h5ad

Example Code

Follow this Jupyter notebook for RAPIDS analysis of this dataset. In order for the notebook to run, the files rapids_scanpy_funcs.py and utils.py need to be in the same folder as the notebook.

We provide a second notebook with the CPU version of this analysis here.

Acceleration

We report the runtime of these notebooks on various GCP instances below. All runtimes are given in seconds. Acceleration is given in parentheses. Benchmarking was performed on Dec 16, 2020. Note: this section is out of date and will be revised shortly.

Step	CPU n1-highmem-64 64 vCPUs	GPU n1-highmem-16 T4 16 GB GPU (Acceleration)	GPU n1-highmem-16 Tesla V100 16 GB GPU (Acceleration)	GPU a2-highgpu-1g Tesla A100 40GB GPU (Acceleration)
Data load + Preprocessing	1120	1125 (1x)	967 (1.2x)	475 (2.4x)
PCA	44	45 (1x)	43 (1x)	17.8 (2.5x)
t-SNE	6509	196 (33x)	50 (130x)	37 (176x)
k-means (single iteration)	148	12.7 (12x)	2.6 (57x)	2 (74x)
KNN	154	141 (1.1x)	92 (1.7x)	62 (2.5x)
UMAP	2571	146 (18x)	32 (80x)	21 (122x)
Louvain clustering	1153	6.1 (189x)	3.9 (296x)	2.4 (480x)
Leiden clustering	6345	5.1 (1244x)	2.7 (2350x)	1.7 (3732x)
Re-analysis of subgroup	255	19.2 (13x)	15 (17x)	17.9 (14.2x)
End-to-end notebook run	18338	1759	1265	686
Price ($/hr)	3.786	1.296	5.906	3.673
Total cost ($)	19.285	0.633	2.075	0.700

Example 3: GPU-based Interactive Visualization of 70,000 Human Lung Cells (beta version)

We demonstrate how to use RAPIDS, Scanpy and Plotly Dash to create an interactive dashboard where we visualize a single-cell RNA-sequencing dataset. Within the interactive dashboard, we can cluster, visualize, and compare any selected groups of cells.

Installation

Additional dependencies are needed for this example. Follow these instructions for conda installation:

conda env create --name rapidgenomics-viz -f conda/rapidgenomics_cuda11.0.viz.yml
conda activate rapidgenomics-viz
python -m ipykernel install --user --display-name "Python (rapidgenomics-viz)"

After installing the necessary dependencies, you can just run jupyter lab.

Example Dataset

The dataset used here is the same as in example 1.

Example Code

Follow this Jupyter notebook to create the interactive visualization. In order for the notebook to run, the files rapids_scanpy_funcs.py and visualize.py need to be in the same folder as the notebook.

Example 4: Droplet Single-cell ATAC-seq of 60K Bone Marrow Cells

We demonstrate the use of RAPIDS to accelerate the analysis of single-cell ATAC-seq data from 60,495 cells. We start with the peak-cell matrix, then perform peak selection, normalization, dimensionality reduction, clustering, and visualization. We also visualize regulatory activity at marker genes and compute differential peaks.

Example Dataset

The dataset is taken from Lareau et al., Nat Biotech 2019. We processed the dataset to include only cells in the 'Resting' condition and peaks with nonzero coverage. Use the following command to download (1) the processed peak-cell count matrix for this dataset (.h5ad), (2) the set of nonzero peak names (.npy), and (3) the cell metadata (.csv), and store them in the data folder:

wget -P <path to this repository>/data https://rapids-single-cell-examples.s3.us-east-2.amazonaws.com/dsci_resting_nonzeropeaks.h5ad; \
wget -P <path to this repository>/data https://rapids-single-cell-examples.s3.us-east-2.amazonaws.com/dsci_resting_peaknames_nonzero.npy; \
wget -P <path to this repository>/data https://rapids-single-cell-examples.s3.us-east-2.amazonaws.com/dsci_resting_cell_metadata.csv

Example Code

Follow this Jupyter notebook for RAPIDS analysis of this dataset. In order for the notebook to run, the files rapids_scanpy_funcs.py and utils.py need to be in the same folder as the notebook.

We provide a second notebook with the CPU version of this analysis here.

Acceleration

We report the runtime of these notebooks on various GCP instances below. All runtimes are given in seconds. Acceleration is given in parentheses. Benchmarking was performed on Dec 16, 2020.

Step	CPU n1-standard-16 16 vCPUs	GPU n1-standard-16 T4 16 GB GPU (Acceleration)	GPU n1-highmem-8 Tesla V100 16 GB GPU (Acceleration)	GPU a2-highgpu-1g Tesla A100 40GB GPU (Acceleration)
PCA	149	146 (1x)	68 (2.2x)	54 (2.8x)
KNN	19.7	19.3 (1x)	5.8 (3.4x)	5.3 (3.7x)
UMAP	69	1.1 (63x)	0.71 (97x)	0.69 (100x)
Louvain clustering	13.1	0.13 (100x)	0.11 (119x)	0.11 (119x)
Leiden clustering	15.7	0.08 (196x)	0.07 (224x)	0.06 (262x)
t-SNE	258	3.2 (81x)	1.5 (172x)	2.2 (117x)
Differential Peak Analysis	135	59 (2.3x)	14.8 (9x)	10.4 (13x)
End-to-end notebook run	682	263	107	92
Price ($/hr)	0.760	1.110	2.953	3.673
Total cost ($)	0.144	0.081	0.110	0.094

Example 5: Visualizing Chromatin Accessibility in 5,000 PBMCs with RAPIDS and AtacWorks (Beta version)

We analyze single-cell ATAC-seq data from 5000 PBMC cells as in example 4. Additionally, we use cuDF to calculate and visualize cluster-specific chromatin accessibility in selected marker regions. Finally, we use a deep learning model trained with AtacWorks, to improve the accuracy of the chromatin accessibility track and call peaks in individual clusters.

Example Data

The dataset was made publicly available by 10X Genomics. Use the following command to download the peak x cell count matrix and the fragment file for this dataset, and store both in the data folder:

wget -P <path to this repository>/data https://rapids-single-cell-examples.s3.us-east-2.amazonaws.com/5k_pbmcs_10X.sparse.h5ad
wget -P <path to this repository>/data https://rapids-single-cell-examples.s3.us-east-2.amazonaws.com/5k_pbmcs_10X_fragments.tsv.gz
wget -P <path to this repository>/data https://rapids-single-cell-examples.s3.us-east-2.amazonaws.com/5k_pbmcs_10X_fragments.tsv.gz.tbi

Example Model

We use a pre-trained deep learning model to denoise the chromatin accessibility track and call peaks. This model can be downloaded into the models folder:

wget -P <path to this repository>/models https://api.ngc.nvidia.com/v2/models/nvidia/atac_bulk_lowcov_5m_50m/versions/0.3/files/models/model.pth.tar

Example Code

Follow this Jupyter notebook for GPU analysis of this dataset. In order for the notebook to run, the files utils.py, and coverage.py need to be in the same folder as the notebook.

Adapting these examples to another dataset

For our examples, we stored the count matrix in a sparse .h5ad format. To convert a different count matrix into this format, follow the instructions in this notebook.