Over the last decade, the size of scRNA-seq datasets has exploded to millions of cells sequenced in a single study. This allows us profile the transcriptomes of various cell types across tissues. The rapid growth of scRNA-seq data has also created an unique set of challenges, for instance, there is a pressing need for scalable approaches for scRNA-seq data visualization.
This vignette introduces scBubbletree1 https://doi.org/10.1186/s12859-024-05927-y, a transparent workflow for quantitative exploration of single cell RNA-seq data.
In short, the algorithm of scBubbletree performs clustering to identify clusters (“bubbles”) of transcriptionally similar cells, and then visualizes these clusters as leafs in a hierarchical dendrogram (“bubbletree”) which describes their relationships. The workflow comprises four steps: 1. determining the clustering resolution, 2. clustering, 3. hierarchical cluster grouping and 4. visualization. We explain each step in the following using scRNA-seq dataset of five cancer cell lines.
To run this vignette we need to load a few packages:
library(scBubbletree)
library(ggplot2)
library(ggtree)
library(patchwork)
Here we will analyze a scRNA-seq data3 https://doi.org/10.1038/s41592-019-0425-8 containing a mixture of 3,918 cells from five human lung adenocarcinoma cell lines (HCC827, H1975, A549, H838 and H2228). The dataset is available here4 https://github.com/LuyiTian/sc_mixology/blob/master/data/ sincell_with_class_5cl.RData.
The library has been prepared with 10x Chromium platform and sequenced with Illumina NextSeq 500 platform. Raw data has been processed with cellranger. The tool demuxlet has been used to predict the identity of each cell based on known genetic differences between the different cell lines.
Data processing was performed with the R-package Seurat. Gene expressions were normalized with the function SCTransform using default parameters, and principal component analysis (PCA) was performed with function RunPCA based on the 5,000 most variable genes in the dataset identified with the function FindVariableFeatures.
In the data we saw that the first 15 principal components capture most of the variance in the data, and the proportion of variance explained by each subsequent principal component was negligible. Thus, we used the single cell projections (embeddings) in 15-dimensional feature space, \(A^{3,918\times 15}\).
# # This script can be used to generate data("d_ccl", package = "scBubbletree")
#
# # create directory
# dir.create(path = "case_study/")
#
# # download the data from:
# https://github.com/LuyiTian/sc_mixology/raw/master/data/
# sincell_with_class_5cl.RData
#
# # load the data
# load(file = "case_study/sincell_with_class_5cl.RData")
#
# # we are only interested in the 10x data object 'sce_sc_10x_5cl_qc'
# d <- sce_sc_10x_5cl_qc
#
# # remove the remaining objects (cleanup)
# rm(sc_Celseq2_5cl_p1, sc_Celseq2_5cl_p2, sc_Celseq2_5cl_p3, sce_sc_10x_5cl_qc)
#
# # get the meta data for each cell
# meta <- colData(d)[,c("cell_line_demuxlet","non_mt_percent","total_features")]
#
# # create Seurat object from the raw counts and append the meta data to it
# d <- Seurat::CreateSeuratObject(counts = d@assays$data$counts,
# project = '')
#
# # check if all cells are matched between d and meta
# # table(rownames(d@meta.data) == meta@rownames)
# d@meta.data <- cbind(d@meta.data, meta@listData)
#
# # cell type predictions are provided as part of the meta data
# table(d@meta.data$cell_line)
#
# # select 5,000 most variable genes
# d <- Seurat::FindVariableFeatures(object = d,
# selection.method = "vst",
# nfeatures = 5000)
#
# # Preprocessing with Seurat: SCT transformation + PCA
# d <- SCTransform(object = d,
# variable.features.n = 5000)
# d <- RunPCA(object = d,
# npcs = 50,
# features = VariableFeatures(object = d))
#
# # perform UMAP + t-SNE
# d <- RunUMAP(d, dims = 1:15)
# d <- RunTSNE(d, dims = 1:15)
#
# # save the preprocessed data
# save(d, file = "case_study/d.RData")
#
# # save the PCA matrix 'A', meta data 'm' and
# # marker genes matrix 'e'
# d <- get(load(file ="case_study/d.RData"))
# A <- d@reductions$pca@cell.embeddings[, 1:15]
# m <- d@meta.data
# e <- t(as.matrix(d@assays$SCT@data[
# rownames(d@assays$SCT@data) %in%
# c("ALDH1A1",
# "PIP4K2C",
# "SLPI",
# "CT45A2",
# "CD74"), ]))
#
# d_ccl <- list(A = A, m = m, e = e)
# save(d_ccl, file = "data/d_ccl.RData")
Load the processed PCA matrix and the meta data
# Load the data
data("d_ccl", package = "scBubbletree")
# Extract the 15-dimensional PCA matrix A
# A has n=cells as rows, f=15 features as columns (e.g. from PCA)
A <- d_ccl$A
dim(A)
#> [1] 3918 15
# Extract the meta-data. For each cell this data contains some
# additional information. Inspect this data now!
m <- d_ccl$m
colnames(m)
#> [1] "orig.ident" "nCount_RNA" "nFeature_RNA"
#> [4] "cell_line_demuxlet" "non_mt_percent" "total_features"
#> [7] "nCount_SCT" "nFeature_SCT"
# Extract the normalized expressions of five marker genes. Rows
# are cells.
e <- d_ccl$e
colnames(e)
#> [1] "ALDH1A1" "SLPI" "CD74" "PIP4K2C" "CT45A2"
We will analyze this data with scBubbletree.
As main input scBubbletree uses matrix \(A^{n\times f}\) which represents a low-dimensional projection of the original scRNA-seq data, with \(n\) rows as cells and \(f\) columns as low-dimension features. Here we use \(A^{3,918\times 15}\) as input.
Important remark about \(A\): the scBubbletree workflow works directly with the numeric matrix \(A^{n\times f}\) and is agnostic to the initial data processing protocol. This enables seamless integration of scBubbletree with computational pipelines using objects generated by the R-packages Seurat and SingleCellExperiment. The users simply have to extract \(A\) from the corresponding Seurat or SingleCellExperiment objects.
The scBubbletree workflow performs the following steps:
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.
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 \(r\)s 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 \approx 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 \(r\)s, 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 \(r\)s 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)
gap_mean | r | k | gap_SE | L95 | H95 |
---|---|---|---|---|---|
2.1686 | 0.0032 | 5 | 0.0064 | 2.1560 | 2.1811 |
2.1639 | 0.0100 | 5 | 0.0039 | 2.1563 | 2.1716 |
2.1678 | 0.0316 | 5 | 0.0048 | 2.1584 | 2.1771 |
2.1653 | 0.1000 | 5 | 0.0084 | 2.1489 | 2.1818 |
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")
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_{\text{eff}}\) (default \(N_{\text{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^{k'\times k'}_{b}\)), 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 (\(\hat{D}^{k' \times k'}\)) and from this a corresponding consensus hierarchical dendrogram (bubbletree; \(\hat{H}\)) 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 (\(\approx\) 32% of all cells in the dataset). Bubble 4 is the smallest one and contains only 437 cells (\(\approx\) 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)
label | Cells | n | p | pct | lab_short | lab_long | tree_order |
---|---|---|---|---|---|---|---|
4 | 437 | 3918 | 0.11 | 11.2 | 4 (0.4K, 11.2%) | 4 (437, 11.2%) | 5 |
3 | 590 | 3918 | 0.15 | 15.1 | 3 (0.6K, 15.1%) | 3 (590, 15.1%) | 4 |
0 | 1253 | 3918 | 0.32 | 32.0 | 0 (1.3K, 32%) | 0 (1253, 32%) | 3 |
2 | 761 | 3918 | 0.19 | 19.4 | 2 (0.8K, 19.4%) | 2 (761, 19.4%) | 2 |
1 | 877 | 3918 | 0.22 | 22.4 | 1 (0.9K, 22.4%) | 1 (877, 22.4%) | 1 |
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.
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")
The two dendrograms are shown side-by-side.
l$tree|k$tree
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_D\)s) 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