Determine the clustering resolution (step 1)

How many clusters (cell types) are there are in the data?

Before we apply clustering, we must first find appropriate value for the resolution parameter k (if we intend to use k-means) or r (if we intend to use graph based community detection approaches such as Louvain). In the next we will first perform Louvain clustering and then k-means clustering.

Determining the resolution parameter r for Louvain clustering

How many clusters (cell types) are there are in the data?

For Louvain clustering we need to select a clustering resolution r. Higher resolutions lead to more communities (k′) and lower resolutions lead to fewer communities.

To find a reasonable value of r we can study the literature or databases such as the human protein atlas database (HPA). We can also use the function get_r for data-driven inference of r based on the Gap statistic.

Lets use the function get_r for data-driven estimation of r based on the Gap statistic and WCSS. As input we need to provide the matrix A and a vector of rs to inspect. See the help function ?get_r to learn more about the remaining input parameters. The output of get_r is the Gap statistic and WCSS estimate for each r (or the number of communities k′ detected at resolution r).

Lets run get_r now [this might take a minute]:

b_r <- get_r(B_gap = 5,
             rs = 10^seq(from = -4, to = 0.5, by = 0.5),
             x = A,
             n_start = 10,
             iter_max = 50,
             algorithm = "original",
             knn_k = 50,
             cores = 1)

The Gap curve has noticeable knee (elbow) at r ≈ 0.003 (dashed gray line). Means (points) and 95% confidence intervals are shown for the Gap statistic at each r using B_gap=5 MCMC simulations.

ggplot(data = b_r$gap_stats_summary)+
  geom_line(aes(x = r, y = gap_mean))+
  geom_point(aes(x = r, y = gap_mean), size = 1)+
  geom_errorbar(aes(x = r, y = gap_mean, ymin = L95, ymax = H95), width = 0.1)+
  ylab(label = "Gap")+
  xlab(label = "r")+
  geom_vline(xintercept = 0.003, col = "gray", linetype = "dashed")+
  scale_x_log10()+
  annotation_logticks(base = 10, sides = "b")

The resolutions r are difficult to interpret. Lets map r to the number of detected communities k′ (analogous to clusters k in k-means clustering), and show the Gap curve as a function of k′.

ggplot(data = b_r$gap_stats_summary)+
  geom_line(aes(x = k, y = gap_mean))+
  geom_point(aes(x = k, y = gap_mean), size = 1)+
  geom_errorbar(aes(x = k, y = gap_mean, ymin = L95, ymax = H95), width = 0.1)+
  geom_vline(xintercept = 5, col = "gray", linetype = "dashed")+
  ylab(label = "Gap")+
  xlab(label = "k'")

A range of resolutions yield k′=5 number of communities, i.e. among the tested rs, we saw k′=5 communities for r = 0.003, 0.01, 0.03 and 0.1. We can e.g. use r=0.1 for our clustering.

ggplot(data = b_r$gap_stats_summary)+
  geom_point(aes(x = r, y = k), size = 1)+
  xlab(label = "r")+
  ylab(label = "k'")+
  scale_x_log10()+
  annotation_logticks(base = 10, sides = "b")+
  theme_bw()

Table with rs that match to k′=5:

knitr::kable(x = b_r$gap_stats_summary[b_r$gap_stats_summary$k == 5, ],
             digits = 4, row.names = FALSE)

Determining the number of clusters k for k-means clustering [alternative]

If we want to use k-means for clustering, then we need to find a reasonable value of k, e.g. by applying once again a data-driven search for k using get_k.

Here get_k will inspect the Gap and WCSS at k = 1, 2, …, 10.

b_k <- get_k(B_gap = 5,
             ks = 1:10, 
             x = A,
             n_start = 50,
             iter_max = 200,
             kmeans_algorithm = "MacQueen",
             cores = 1)

Notice the similar Gap curve with noticeable knee (elbow) at k = 5 (dashed gray line). Means (points) and 95% confidence intervals are shown for the Gap statistic at each k using B_gap=5 MCMC simulations.

ggplot(data = b_k$gap_stats_summary)+
  geom_line(aes(x = k, y = gap_mean))+
  geom_point(aes(x = k, y = gap_mean), size = 1)+
  geom_errorbar(aes(x = k, y = gap_mean, ymin = L95, ymax = H95), width = 0.1)+
  ylab(label = "Gap")+
  geom_vline(xintercept = 5, col = "gray", linetype = "dashed")

Clustering with Louvain

Now that we found out that r = 0.1 (k′ = 5) is a reasonable choice based on the data, we will perform Louvain clustering with r = 0.1 and A as inputs. For this we will use the function get_bubbletree_graph.

After the clustering is complete we will organize the bubbles by hierarchical clustering. For this we perform B bootstrap iterations. In iteration b the algorithm draws a random subset of N_eff (default N_eff = 200) cells with replacement from each cluster and compute the average inter-cluster Euclidean distances. This data is used to populate the distance matrix (D_b^{k′ × k′}), which is provided as input for hierarchical clustering with average linkage to generate a hierarchical clustering dendrogram H_b.

The collection of distance matrices that are computed during B iterations are used to compute a consensus (average) distance matrix (D̂^{k′ × k′}) and from this a corresponding consensus hierarchical dendrogram (bubbletree; Ĥ) is constructed. The collection of dendrograms are used to quantify the robustness of the bubbletree topology, i.e. to count the number of times each branch in the bubbletree is found among the topologies of the bootstrap dendrograms. Branches can have has variable degrees of support ranging between 0 (no support) and B (complete support). Distances between bubbles (inter- bubble relationships) are described quantitatively in the bubbletree as sums of branch lengths.

Steps 2. (clustering) and 3. (hierarchical grouping) are performed now:

l <- get_bubbletree_graph(x = A,
                          r = 0.1,
                          algorithm = "original",
                          n_start = 20,
                          iter_max = 100,
                          knn_k = 50,
                          cores = 1,
                          B = 300,
                          N_eff = 200,
                          round_digits = 1,
                          show_simple_count = FALSE)

#  See the help `?get_bubbletree_graph` to learn about the input parameters.

… and plot the bubbletree

l$tree

Lets describe the bubbletree:

bubbles: The bubbletree has k'=5 bubbles (clusters) shown as leaves. The absolute and relative cell frequencies in each bubble and the bubble IDs are shown as labels. Bubble radii scale linearly with absolute cell count in each bubble, i.e. large bubbles have many cells and small bubbles contain few cells.

Bubble 0 is the largest one in the dendrogram and contains 1,253 cells (≈ 32% of all cells in the dataset). Bubble 4 is the smallest one and contains only 437 cells (≈ 11% of all cells in the dataset).

We can access the bubble data shown in the bubbletree

knitr::kable(l$tree_meta, digits = 2, row.names = FALSE)

topology: inter-bubble distances are represented by sums of branch lengths in the dendrogram. Branches of the bubbletree are annotated with their bootstrap support values (red branch labels). The branch support value tells us how many times a given branch from the bubbletree was found among the B bootstrap dendrograms. We ran get_bubbletree_graph with B = 300. All but one branch have complete (300 out of 300) support, and one branch has lower support of 270 (90%). This tells us that the branch between bubbles (3, 4) and 0 is not as robust.

Clustering with k-means ([alternative])

To perform clustering with the k-means method we can use the function get_bubbletree_kmeans.

k <- get_bubbletree_kmeans(x = A,
                           k = 5,
                           cores = 1,
                           B = 300,
                           N_eff = 200,
                           round_digits = 1,
                           show_simple_count = FALSE,
                           kmeans_algorithm = "MacQueen")

Comparison of bubbletree based on Louvain and k-means clustering

The two dendrograms are shown side-by-side.

l$tree|k$tree

Visual comparison of two bubbletrees

To compare a pair of bubbletrees generated based on the same data but with different inputs we can use the function compare_bubbletrees.

The function generates two bubbletrees and a heatmap, where the tiles of the heatmap are color coded according to the jaccard distances (J_Ds) between the pairs of bubbles from the two bubbletrees, and the tile labels show the numbers of cells in common between the bubbles.

Reminder of the Jaccard index (J) and the Jaccard distance (J_D):

For clusters A and B we compute the Jaccard index $J(A,B)=\dfrac{|A \cap B|}{|A \cup B|}$ and distance J_D(A, B) = 1 − J(A, B). If A and B contain the same set of cells J_D(A, B)=0, and if they have no cells in common J_D(A, B)=1.

The heatmap hints at nearly identical clusterings between the bubbletrees. Only 3 cells (red tiles with label = 1) are classified differently by the two bubbletrees.

cp <- compare_bubbletrees(btd_1 = l, 
                          btd_2 = k,
                          ratio_heatmap = 0.6,
                          tile_bw = F,
                          tile_text_size = 3)
cp$comparison

Attaching categorical features

Categorical cell features can be ‘attached’ to the bubbletree using the function get_cat_tiles. Here we will show the relative frequency of cell type labels across the bubbles (parameter integrate_vertical=TRUE).

Interpretation of the figure below:

• we see high degree of co-occurrence between cell lines and bubbles, i.e.  each bubble is made up of cells from a distinct cell line
• for instance, 99.8% of cells that have feature HCC827 are found in bubble 3
• columns in the tile plot integrate to 100%

w1 <- get_cat_tiles(btd = l,
                    f = m$cell_line_demuxlet,
                    integrate_vertical = TRUE,
                    round_digits = 1,
                    x_axis_name = 'Cell line',
                    rotate_x_axis_labels = TRUE,
                    tile_text_size = 2.75)

(l$tree|w1$plot)+
  patchwork::plot_layout(widths = c(1, 1))

We can also show the inter-bubble cell type composition, i.e. the relative frequencies of different cell types in a specific bubble (with parameter integrate_vertical=FALSE).

Interpretation of the figure below:

• the bubbles appear to be “pure” → made up of cells from distinct cell lines
• the cell line composition of bubble 2 is: 0.1% H838, 99.6% H2228, 0.1% A549, 0.1% H1975 and 0% HCC827 cells
• rows integrate to 100% (here the numbers in the heatmap tiles are rounded to the nearest tenth, hence they integrate approximately to 100%)

w2 <- get_cat_tiles(btd = l,
                    f = m$cell_line_demuxlet,
                    integrate_vertical = FALSE,
                    round_digits = 1,
                    x_axis_name = 'Cell line',
                    rotate_x_axis_labels = TRUE,
                    tile_text_size = 2.75)

(l$tree|w2$plot)+patchwork::plot_layout(widths = c(1, 1))

scBubbletree uses R-package ggtree to visualize the bubbletree, and ggplot2 to visualize annotations. Furthermore, R-package patchwork is used to combine plots.

(l$tree|w1$plot|w2$plot)+
  patchwork::plot_layout(widths = c(1, 2, 2))+
  patchwork::plot_annotation(tag_levels = "A")

Gini impurity index

To quantify the purity of a cluster (or bubble) i with n_i number of cells, each of which carries one of L possible labels (e.g. cell lines), we can compute the Gini impurity index:

$\textit{GI}_i=\sum_{j=1}^{L} \pi_{ij}(1-\pi_{ij})$,

with π_ij as the relative frequency of label j in cluster i. In homogeneous (pure) clusters most cells carry a distinct label. Hence, the π’s are close to either 1 or 0, and GI takes on a small value close to zero. In impure clusters cells carry a mixture of different labels. In this case most π are far from either 1 or 0, and GI diverges from 0 and approaches 1. If the relative frequencies of the different labels in cluster i are equal to the (background) relative frequencies of the labels in the sample, then cluster i is completely impure.

To compute the overall Gini impurity of a bubbletree, which represents a clustering solution with k bubbles, we estimated the weighted Gini impurity (WGI) by computing the weighted (by the cluster size) average of the

$\textit{WGI}=\sum_{i=1}^{k} \textit{GI}_i \dfrac{n_i}{n}$,

with n_i as the number of cells in cluster i and n = ∑_in_i.

The Gini impurity results are shown below:

# gini
get_gini(labels = m$cell_line_demuxlet, 
         clusters = l$cluster)$gi

#>   cluster          GI
#> 1       3 0.010129273
#> 2       1 0.004553202
#> 3       4 0.000000000
#> 4       2 0.007863642
#> 5       0 0.000000000

All cluster-specific GIs are close to 0 and hence also WGI is close to 0. This indicates nearly perfect mapping of cell lines to bubbles. This analysis performed for different values of r with function get_gini_r, which takes as main input the output of get_k or get_r

gini_boot <- get_gini_k(labels = m$cell_line_demuxlet, obj = b_r)

From the figure we can conclude that WGI drops to 0 at k=5, and all labels are nearly perfectly split across the bubbles with each bubble containing cells exclusively from one cell type.

g1 <- ggplot(data = gini_boot$wgi_summary)+
  geom_point(aes(x = k, y = wgi), size = 1)+
  ylab(label = "WGI")+
  ylim(c(0, 1))

g1

Attaching numeric features

We can also “attach” numeric cell features to the bubbletree. We will “attach” the expression of five marker genes, i.e. one marker gene for each of the five cancer cell lines.

We can visualize numeric features in two ways.

First, we can show numeric feature aggregates (e.g. “mean”, “median”, “sum”, “pct nonzero” or “pct zero”) in the different bubbles with get_num_tiles

w3 <- get_num_tiles(btd = l,
                    fs = e,
                    summary_function = "mean",
                    x_axis_name = 'Gene expression',
                    rotate_x_axis_labels = TRUE,
                    round_digits = 1,
                    tile_text_size = 2.75)

(l$tree|w3$plot)+patchwork::plot_layout(widths = c(1, 1))

Second, we can visualize the distributions of the numeric cell features in each bubble as violins with get_num_violins

w4 <- get_num_violins(btd = l,
                      fs = e,
                      x_axis_name = 'Gene expression',
                      rotate_x_axis_labels = TRUE)

(l$tree|w3$plot|w4$plot)+
  patchwork::plot_layout(widths = c(1.5, 2, 2.5))+
  patchwork::plot_annotation(tag_levels = 'A')

Quality control with scBubbletree

What is the percent of UMIs coming from mitochondrial genes in each bubble?

w_mt_dist <- get_num_violins(btd = l,
                             fs = 1-m$non_mt_percent,
                             x_axis_name = 'MT [%]',
                             rotate_x_axis_labels = TRUE)

w_umi_dist <- get_num_violins(btd = l,
                              fs = m$nCount_RNA/1000,
                              x_axis_name = 'RNA count (in thousands)',
                              rotate_x_axis_labels = TRUE)

w_gene_dist <- get_num_violins(btd = l,
                               fs = m$nFeature_RNA,
                               x_axis_name = 'Gene count',
                               rotate_x_axis_labels = TRUE)


(l$tree|w_mt_dist$plot|w_umi_dist$plot|w_gene_dist$plot)+
  patchwork::plot_layout(widths = c(1, 1, 1, 1))+
  patchwork::plot_annotation(tag_levels = 'A')