library(knitr)
opts_chunk$set(fig.align = 'center', fig.width = 6, fig.height = 5)

Introduction

ILoReg is a tool for cell population identification from single-cell RNA-seq (scRNA-seq) data. In our paper [1], we showed that ILoReg was able to identify, by both unsupervised clustering and visually, rare cell populations that other scRNA-seq data analysis pipelines were unable to identify.

The figure below illustrates the workflows of ILoReg and a typical pipeline that applies feature selection prior to dimensionality reduction by principal component analysis (PCA).

*Figure: Workflows of ILoReg and a feature-selection based approach.*{width=100%}

In contrast to most scRNA-seq data analysis pipelines, ILoReg does not reduce the dimensionality of the gene expression matrix by feature selection. Instead, it performs probabilistic feature extraction using iterative clustering projection (ICP), generating an $N \times k$ -dimensional probability matrix, which contains probabilities of each of the $N$ cells belonging to the $k$ clusters. ICP is a novel self-supervised learning algorithm that iteratively seeks a clustering with $k$ clusters that maximizes the adjusted Rand index (ARI) between the clustering $C$ and its projection $C'$ by L1-regularized logistic regression. In the ILoReg consensus approach, ICP is run $L$ times and the $L$ probability matrices are merged into a joint probability matrix and subsequently transformed by principal component analysis (PCA) into a lower dimensional ($N \times p$) matrix (consensus matrix). The final clustering step is performed using hierarhical clustering by the Ward's method, after which the user can extract a clustering with $K$ consensus clusters. However, the user can also use any other clustering method at this point. Two-dimensional visualization is supported using two popular nonlinear dimensionality reduction methods: t-distributed stochastic neighbor embedding (t-SNE) and uniform manifold approximation and projection (UMAP). Additionally, ILoReg provides user-friendly functions that enable identification of differentially expressed (DE) genes and visualization of gene expression.

Installation

ILoReg can be downloaded from Bioconductor and installed by executing the following command in the R console.

if(!requireNamespace("BiocManager", quietly = TRUE))
    install.packages("BiocManager")
BiocManager::install("ILoReg")

Example: Peripheral Blood Mononuclear Cells

In the following, we go through the different steps of ILoReg's workflow and demonstrate it using a peripheral blood mononuclear cell (PBMC) dataset. The toy dataset included in the ILoReg R package (pbmc3k_500) contains 500 cells that have been downsampled from the pbmc3k dataset [2]. The preprocessing was rerun with a newer reference genome (GRCh38.p12) and Cell Ranger v2.2.0 [3] to identify different immunoglobulin subpopulations in B-cells.

Setup a SingleCellExperiment object and prepare it for ILoReg analysis

The only required input for ILoReg is a log-transformed, normalized gene expression matrix that has been, with genes/features in rows and cells/samples in columns. The input can be of matrix, data.frame or dgCMatrix class, which is then transformed into a sparse object of dgCMatrix class. Please note that the method has been designed to work with sparse data, i.e. with a high proportion of zero values. If, for example, the features of your dataset have been standardized, the run time and the memory usage of ILoReg will likely be much higher.

suppressMessages(library(ILoReg))
suppressMessages(library(SingleCellExperiment))
suppressMessages(library(cowplot))
# The dataset was normalized using the LogNormalize method from the Seurat R package.
sce <- SingleCellExperiment(assays = list(logcounts = pbmc3k_500))
sce <- PrepareILoReg(sce)

Run the ICP clustering algorithm $L$ times

Running ICP $L$ times in parallel is the most computationally demanding part of the workflow.

In the following, we give a brief summary of the parameters.

As general guidelines on how to adjust the parameters, we recommend leaving $r$ and $L$ to their defaults ($r=5$ and $L=200$). However, increasing $k$ from 15 to e.g. 30 can reveal new cell subsets that are of interest. Regarding $d$, increasing it to somewhere between 0.4-0.6 helps if the user wants lower resolution (less distinguishable populations). Increasing $C$ from 0.3 to 1 reduces the number of distinguishable populations, as the logistic regression model filters out fewer genes.

# ICP is stochastic. To obtain reproducible results, use set.seed().
set.seed(1)
# Run ICP L times. This is  the slowest step of the workflow, 
# and parallel processing can be used to greatly speed it up.
sce <- RunParallelICP(object = sce, k = 15,
                      d = 0.3, L = 30, 
                      r = 5, C = 0.3, 
                      reg.type = "L1", threads = 0)

PCA transformation of the joint probability matrix

The $L$ probability matrices are merged into a joint probability matrix, which is then transformed into a lower dimensionality by PCA. Before applying PCA, the user can optionally scale the cluster probabilities to unit-variance.

# p = number of principal components
sce <- RunPCA(sce,p=50,scale = FALSE)

Optional: PCA requires the user to specify the number of principal components, for which we selected the default value $p=50$. To aid in decision making, the elbow plot is commonly used to seek an elbow point, of which proximity the user selects $p$. In this case the point would be close to $p=10$. Trying both a $p$ that is close to the elbow point and the default $p=50$ is recommended.

PCAElbowPlot(sce)

Nonlinear dimensionality reduction

To visualize the data in two-dimensional space, nonlinear dimensionality reduction is performed using t-SNE or UMAP. The input data for this step is the $N \times p$ -dimensional consensus matrix.

sce <- RunUMAP(sce)
sce <- RunTSNE(sce,perplexity=30)

Gene expression visualization

Visualize the t-SNE and UMAP transformations using the GeneScatterPlot function, highlighting expression levels of CD3D (T cells), CD79A (B cells), CST3 (monocytes, dendritic cells, platelets), FCER1A (myeloid dendritic cells).

GeneScatterPlot(sce,c("CD3D","CD79A","CST3","FCER1A"),
                dim.reduction.type = "umap",
                point.size = 0.3)
GeneScatterPlot(sce,c("CD3D","CD79A","CST3","FCER1A"),
                dim.reduction.type = "tsne",
                point.size = 0.3)

Hierarchical clustering using the Ward's method

The $N \times p$ -dimensional consensus matrix is hierarchically clustered using the Ward's method.

sce <- HierarchicalClustering(sce)

Extracting a consensus clustering with $K$ clusters

After the hierarchical clustering, the user needs to define how many consensus clusters ($K$) to extract from the tree dendrogram. The SelectKClusters function enables extracting a consensus clustering with $K$ clusters. Please note that the clustering is overwritten every time the function is called.

# Extract K=13 clusters.
sce <- SelectKClusters(sce,K=13)

Next, we use the ClusteringScatterPlot function to draw the t-SNE and UMAP transformations and color each cell according to the cluster labels.

# Use plot_grid function from the cowplot R package to combine the two plots into one.
plot_grid(ClusteringScatterPlot(sce,
                                dim.reduction.type = "umap",
                                return.plot = TRUE,
                                title = "UMAP",
                                show.legend=FALSE),
          ClusteringScatterPlot(sce,
                                dim.reduction.type = "tsne",
                                return.plot = TRUE
                                ,title="t-SNE",
                                show.legend=FALSE),
          ncol = 1
)

Identification of gene markers

TheILoReg R package provides functions for the identification of gene markers of clusters. This is accomplished by DE analysis, where gene expression levels of the cells from each cluster are compared against the rest of the cells. Currently, the only supported statistical test is the the Wilcoxon rank-sum test (aka Mann-Whitney U test). The p-values are corrected for multiple comparisons using the Bonferroni method. To accelerate the analysis, genes that are less likely to be DE can be filtered out prior to the statistical testing using multiple criteria.

gene_markers <- FindAllGeneMarkers(sce,
                                   clustering.type = "manual",
                                   test = "wilcox",
                                   log2fc.threshold = 0.25,
                                   min.pct = 0.25,
                                   min.diff.pct = NULL,
                                   min.cells.group = 3,
                                   return.thresh = 0.01,
                                   only.pos = TRUE,
                                   max.cells.per.cluster = NULL)

Selecting top gene markers

Select top 10 and 1 gene markers based on the log2 fold-change and the Bonferroni adjusted p-value.

top10_log2FC <- SelectTopGenes(gene_markers,
                               top.N = 10,
                               criterion.type = "log2FC",
                               inverse = FALSE)
top1_log2FC <- SelectTopGenes(gene_markers,
                              top.N = 1,
                              criterion.type = "log2FC",
                              inverse = FALSE)
top10_adj.p.value <- SelectTopGenes(gene_markers,
                                    top.N = 10,
                                    criterion.type = "adj.p.value",
                                    inverse = TRUE)
top1_adj.p.value <- SelectTopGenes(gene_markers,
                                   top.N = 1,
                                   criterion.type = "adj.p.value",
                                   inverse = TRUE)

Draw the t-SNE and UMAP transformations, highlighting expression levels of the top 1 gene markers based on the log2 fold-change.

GeneScatterPlot(sce,
                genes = unique(top1_log2FC$gene),
                dim.reduction.type = "tsne",
                point.size = 0.5,ncol=2)

Gene marker heatmap

GeneHeatmap function enables visualizing gene markers in a heatmap, where cells and genes are grouped by the clustering.

GeneHeatmap(sce,
            clustering.type = "manual",
            gene.markers = top10_log2FC)

Renaming clusters

RenameAllClusters enables renaming all clusters at once.

sce <- RenameAllClusters(sce,
                         new.cluster.names = c("GZMK+/CD8+ T cells",
                                               "IGKC+ B cells",
                                               "Naive CD4+ T cells",
                                               "NK cells",
                                               "CD16+ monocytes",
                                               "CD8+ T cells",
                                               "CD14+ monocytes",
                                               "IGLC+ B cells",
                                               "Intermediate monocytes",
                                               "IGKC+/IGLC+ B cells",
                                               "Memory CD4+ T cells",
                                               "Naive CD8+ T cells",
                                               "Dendritic cells"))

Draw the gene heatmap again, but with the clusters renamed.

GeneHeatmap(sce,gene.markers = top10_log2FC)

Violin plot visualization

VlnPlot enables visualization of gene expression with cells grouped by clustering.

# Visualize CD3D: a marker of T cells
VlnPlot(sce,genes = c("CD3D"),return.plot = FALSE,rotate.x.axis.labels = TRUE)

Additional functionality

ILoReg provides additional functionality for performing tasks, which are sometimes required in scRNA-seq data analysis.

Estimating the optimal number of clusters

The optimal number of clusters can be estimated by calculating the average silhouette value across the cells for a set of clusterings within a range of different $K$ values (e.g. $[2,50]$). The clustering mathing to the highest average silhouette is saved to clustering.optimal slot. Therefore, the clustering acquired using the SelectKClusters function is not overwritten.

sce <- CalcSilhInfo(sce,K.start = 2,K.end = 50)
SilhouetteCurve(sce,return.plot = FALSE)

Renaming one cluster

sce <- SelectKClusters(sce,K=20)
# Rename cluster 1 as A
sce <- RenameCluster(sce,old.cluster.name = 1,new.cluster.name = "A")

Visualize with a custom annotation

# Select a clustering with K=5 clusters
sce <- SelectKClusters(sce,K=5)
# Generate a custom annotation with K=5 clusters and change the names to the five first alphabets.
custom_annotation <- plyr::mapvalues(metadata(sce)$iloreg$clustering.manual,c(1,2,3,4,5),LETTERS[1:5])
# Visualize the annotation
AnnotationScatterPlot(sce,
                      annotation = custom_annotation,
                      return.plot = FALSE,
                      dim.reduction.type = "tsne",
                      show.legend = FALSE)

Merging clusters

# Merge clusters 3 and 4
sce <- SelectKClusters(sce,K=20)
sce <- MergeClusters(sce,clusters.to.merge  = c(3,4))

Identification of differentially expressed genes between two arbitrary sets of clusters

sce <- SelectKClusters(sce,K=13)
sce <- RenameAllClusters(sce,
                         new.cluster.names = c("GZMK+/CD8+ T cells",
                                               "IGKC+ B cells",
                                               "Naive CD4+ T cells",
                                               "NK cells",
                                               "CD16+ monocytes",
                                               "CD8+ T cells",
                                               "CD14+ monocytes",
                                               "IGLC+ B cells",
                                               "Intermediate monocytes",
                                               "IGKC+/IGLC+ B cells",
                                               "Memory CD4+ T cells",
                                               "Naive CD8+ T cells",
                                               "Dendritic cells"))
# Identify DE genes between naive and memory CD4+ T cells
GM_naive_memory_CD4 <- FindGeneMarkers(sce,
                                       clusters.1 = "Naive CD4+ T cells",
                                       clusters.2 = "Memory CD4+ T cells",
                                       logfc.threshold = 0.25,
                                       min.pct = 0.25,
                                       return.thresh = 0.01,
                                       only.pos = TRUE)


# Find gene markers for dendritic cells
GM_dendritic <- FindGeneMarkers(sce,
                                clusters.1 = "Dendritic cells",
                                logfc.threshold = 0.25,
                                min.pct = 0.25,
                                return.thresh = 0.01,
                                only.pos = TRUE)

Session info

sessionInfo()

References

  1. Johannes Smolander, Sini Junttila, Mikko S Venäläinen, Laura L Elo, ILoReg: a tool for high-resolution cell population identification from single-cell RNA-seq data, Bioinformatics, Volume 37, Issue 8, 15 April 2021, Pages 1107–1114, https://doi.org/10.1093/bioinformatics/btaa919
  2. "3k PBMCs from a Healthy Donor". 10X Genomics.
  3. "What is Cell Ranger?" 10X Genomics.

Contact information

If you have questions related to ILoReg, please contact us here.



elolab/ILoReg documentation built on March 28, 2022, 1:17 a.m.