Basic Usage

Example of training an HDBSCAN model using the hdbscan Python package in Scikit-learn contrib:

And the same code using the GPU-Accelerated HDBSCAN in cuML (spoiler alert: the only difference is the import).

Plotting

We can plot the minimum spanning tree the same way we would for the original HDBSCAN implementation:

The single linkage dendrogram and condensed tree can be plotted as well:

Clustering Algorithm Comparison

Example notebook showing the strengths of density-based clustering techniques DBSCAN & HDBSCAN on datasets with odd and interleaved shapes.

Generate the data

K-Means

Even when we know there are only 2 clusters in this dataset, we can see that k-means just separates the space into two evenly-sized regions. It's not possible to model the interleaved and odd cluster shapes without performing some heavy pre-processing to project the points into a space where they can be modeled with fixed-sized balls.

DBSCAN

While this particular example dataset tends to be a great demonstration of how DBSCAN can find clusters of odd shapes and similar densities, theepsparameter sets the radius of a ball around each point where points inside the ball are considered part of its neighborhood. This neighborhood is called an epsilon neighborhoods and it can be pretty non-intuitive to tune in practice. In addition, it's assumed that the sameepssetting applies globally to all the data points. A little bit more intuitive is themin_samplesparameter, which determines the neighborhood size within theepsball for a point to be considered acore point, upon which its neighborhood will be expanded to form a cluster.

As we can see, the default paramter settings of DBSCAN (and other density-based algorithms) make assumptions about the data that won't often lead to good clustering results.

We can set themin_samplesparameter to 10, but the amount of noise is too high for this alone to result in good cluster assignments.

We can introduce a setting for theepsparameter, but it's unclear how to determine a good setting for this value just from this visualization. Setting this value too small can yield too many small clusterings which are not representative of the larger patterns in the data.

Increasingepsto 0.12 did the trick and we notice it found the points which surround the dense regions as noise (these will have a label of -1)

HDBSCAN

HDBSCAN removes the need to specify a global epsilon neighborhood radius around each point (eps) and instead uses themin_pointsargument to create a radius to its k-th nearest neighbor, using the radius to push sparser regions further away from each other. Since HDBSCAN is built upon the agglomerative clustering technique single-linkage clustering, each point starts as its own cluster and is merged into a tree, bottom-up, until a single root cluster is reached. Amin_cluster_sizeargument allows us to specify a minimum threshold for when clusters in the agglomerated tree should be merged. The dendrogram is cut at varying levels, or selected, using a measure of cluster stability, based on the distances between the points in each cluster and the clusters from which they originated. Clusters which are selected cosume all of their descendent clusters. HDBSCAN provides an additional optioncluster_selection_epsilonto set the minimum distance threshold from which clusters will be split up into smaller clusters.

Adding the same clustering against the reference Python HDBSCAN implementation

Clustering and Visualizing Embeddings with HDBSCAN & UMAP

In this notebook, we use a dataset of word embeddings to extract areas that could be associated with meaningful topics. We use HDBSCAN to find the meaningful topics since it can account for noise when labeling regions of high density without explicitly requiring a distance as input. We also use UMAP to visualize the resulting topic clusters.

This notebook requires theGoogleNews Word2Vec dataset, which can be downloaded fromhttps://drive.google.com/file/d/0B7XkCwpI5KDYNlNUTTlSS21pQmM/edit?usp=sharing.

We can usemin_samplesto control the radius of the core distances. A smaller setting leads to more clusters and less points being treated as noise.min_cluster_sizeandmax_cluster_sizebound the number of clusters by only allowing points to form a cluster when they are >min_cluster_clusterand splitting a single cluster into multiple smaller clusters when it is >max_cluster_size

L2 normalize the vectors so that the Euclidean distance between them will give the angular distance.

First, we can use manifold learning to convert the neighborhoods of points in the angular distance space to a 2-dimensional Euclidean space so that we can better virualize the clusters.

Percentage of points were labeled as meaningful topics

How many labels are there?

We observe that the regions of high density where points are collected closely together have filled out nicely into clusters of various different shapes and sizes. We interpret the distances of points within their local neighborhoods in the following plot as being more similar than points in disconnected neighborhoods. It's not uncommon for UMAP to split neighborhoods that into multiple clusters even though HDBSCAN labeled them as a single cluster. This just implies that the neighborhoods of the points which are split are more similar to each other internally than they are across their multiple separated components.

We can also visualize the noise

Let's look at a few groups of terms that clustered together

We can also order them according to their probabilities (notice most probabilities are fairly large because the regions are already pretty dense)

RAPIDS & Scanpy Single-Cell RNA-seq Workflow

Copyright (c) 2020, NVIDIA CORPORATION.

Licensed under the Apache License, Version 2.0 (the "License") you may not use this file except in compliance with the License. You may obtain a copy of the License at

http://www.apache.org/licenses/LICENSE-2.0 

Unless required by applicable law or agreed to in writing, software distributed under the License is distributed on an "AS IS" BASIS, WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied. See the License for the specific language governing permissions and limitations under the License.

This notebook demonstrates a single-cell RNA analysis workflow that begins with preprocessing a count matrix of size(n_gene, n_cell)and results in a visualization of the clustered cells for further analysis.

For demonstration purposes, we use a dataset of ~70,000 human lung cells from Travaglini et al. 2020 (https://www.biorxiv.org/content/10.1101/742320v2) and label cells using the ACE2 and TMPRSS2 genes. See the README for instructions to download this dataset.

Import requirements

We use the RAPIDS memory manager on the GPU to control how memory is allocated.

Input data

In the cell below, we provide the path to the.h5adfile containing the count matrix to analyze. Please see the README for instructions on how to download the dataset we use here.

We recommend saving count matrices in the sparse .h5ad format as it is much faster to load than a dense CSV file. To run this notebook using your own dataset, please see the README for instructions to convert your own count matrix into this format. Then, replace the path in the cell below with the path to your generated.h5adfile.

Set parameters

Load and Prepare Data

We load the sparse count matrix from anh5adfile using Scanpy. The sparse count matrix will then be placed on the GPU.

We maintain the index of unique genes in our dataset:

Verify the shape of the resulting sparse matrix:

And the number of non-zero values in the matrix:

Preprocessing

Filter

We filter the count matrix to remove cells with an extreme number of genes expressed.

Some genes will now have zero expression in all cells. We filter out such genes.

The size of our count matrix is now reduced.

Normalize

We normalize the count matrix so that the total counts in each cell sum to 1e4.

Next, we log transform the count matrix.

Select Most Variable Genes

We will now select the most variable genes in the dataset. However, we first save the 'raw' expression values of the ACE2 and TMPRSS2 genes to use for labeling cells afterward. We will also store the expression of an epithelial marker gene (EPCAM).

Now, we convert the count matrix to an annData object.

Using scanpy, we filter the count matrix to retain only the 5000 most variable genes.

Regress out confounding factors (number of counts, ribosomal gene expression)

We can now perform regression on the count matrix to correct for confounding factors - for example purposes, we use the number of counts and the expression of ribosomal genes. Many workflows use the expression of mitochondrial genes (named starting withMT-).

We now calculate the total counts and the percentage of ribosomal counts for each cell.

And perform regression:

Scale

Finally, we scale the count matrix to obtain a z-score and apply a cutoff value of 10 standard deviations, obtaining the preprocessed count matrix.

Cluster & Visualize

We store the preprocessed count matrix as an AnnData object, which is currently in host memory. We also add the expression levels of the marker genes as observations to the annData object.

Reduce

We use PCA to reduce the dimensionality of the matrix to its top 50 principal components.

UMAP + Density clustering

We can also visualize the cells using the UMAP algorithm in Rapids. Before UMAP, we need to construct a k-nearest neighbors graph in which each cell is connected to its nearest neighbors. This can be done conveniently using rapids functionality already integrated into Scanpy.

Note that Scanpy uses an approximation to the nearest neighbors on the CPU while the GPU version performs an exact search. While both methods are known to yield useful results, some differences in the resulting visualization and clusters can be observed.

The UMAP function from Rapids is also integrated into Scanpy.

Next, we use the Louvain algorithm for graph-based clustering, once again using therapidsoption in Scanpy.

We plot the cells using the UMAP visualization, and using the Louvain clusters as labels.