Contents


1 Brief introduction

In this tutorial, we will analyze two datasets: one from Zheng et al., (Nature Communications, 2016) and the other from Biase et al., (Genome Research, 2014). Zheng dataset contains 500 human peripheral blood mononuclear cells (PBMCs) sequenced using GemCode platform, which consists of three cell types, CD56+ natural killer cells, CD19+ B cells and CD4+/CD25+ regulatory T cells. The original data can be downloaded from 10X GENOMICS website. The Biase dataset has 49 mouse embryo cells, which were sequenced by SMART-Seq and can be found at NCBI GEO:GSE57249.


2 Setup the library

library("SAFEclustering")
data("data_SAFE")

3 Zheng dataset

3.1 Setup the input expression matrix

dim(data_SAFE$Zheng.expr)
## [1] 32738   500
data_SAFE$Zheng.expr[1:5, 1:5]
##              CTACAACTCATACG CAACGAACTGGTTG AACGCCCTTTTGCT TATGTGCTAGTGTC
## MIR1302-10                0              0              0              0
## FAM138A                   0              0              0              0
## OR4F5                     0              0              0              0
## RP11-34P13.7              0              0              0              0
## RP11-34P13.8              0              0              0              0
##              CTAAGGTGTTTGCT
## MIR1302-10                0
## FAM138A                   0
## OR4F5                     0
## RP11-34P13.7              0
## RP11-34P13.8              0

3.2 Perform individual clustering

Here we perform single-cell clustering using four popular methods, SC3, CIDR, Seurat and t-SNE + k-means, without filtering any genes or cells.

cluster.results <- individual_clustering(inputTags = data_SAFE$Zheng.expr, datatype = "count", mt_filter = FALSE, nGene_filter = FALSE, SC3 = TRUE, gene_filter = FALSE, CIDR = TRUE, nPC.cidr = NULL, Seurat = TRUE, nPC.seurat = NULL, resolution = 0.9, tSNE = TRUE, dimensions = 3, perplexity = 30, SEED = 123)
## Performing SC3 clustering...
## Estimating k...
## Setting SC3 parameters...
## Calculating distances between the cells...
## Performing transformations and calculating eigenvectors...
## Performing k-means clustering...
## Calculating consensus matrix...
## Performing CIDR clustering...
## Performing Seurat clustering...
## Regressing out: nUMI
## Scaling data matrix
## Performing tSNE + k-means clustering...

The function indiviual_clustering will output a matrix, where each row represents the cluster results of each method, and each colunm represents a cell. User can also extend SAFE-clustering to other scRNA-seq clustering methods, by putting all clustering results into a \(M * N\) matrix with M clustering methods and N cells.

cluster.results[1:4, 1:10]
##      [,1] [,2] [,3] [,4] [,5] [,6] [,7] [,8] [,9] [,10]
## [1,]    2    2    2    2    2    2    2    2    2     2
## [2,]    1    1    1    1    1    1    1    1    1     1
## [3,]    3    3    3    3    3    3    3    3    3     3
## [4,]    1    1    1    1    1    1    1    1    1     1

3.3 Cluster ensemble

Using the clustering results generated in last step, we perform cluster ensemble using three partitioning algorithms meta-clustering algorithm (MCLA), hypergraph partitioning algorithm (HGPA) and cluster-based similarity partitioning algorithm (CSPA) (Strehl and Ghosh, Proceedings of AAAI 2002, Edmonto, Canada, 2002).

Note that HGPA is performed using the shmetis program (from the hMETIS package v. 1.5 (Karypis et al., IEEE Transactions on Very Large Scale Integration (VLSI) Systems, 1999)), and MCLA and CSPA are performed using gpmetis program (from METIS v. 5.1.0 (Karypis and Kumar, SIAM Journal on Scientific Computing, 1998)). Please put them in the working directory or provide the directory where these two programs are.

cluster.ensemble <- SAFE(cluster_results = cluster.results, program.dir = "~/Documents/single_cell_clustering", MCLA = TRUE, CSPA = TRUE, HGPA = TRUE, SEED = 123)

Here is the list of ANMI results for esemble solution of each K and each partitioning algorithm.

## [1] "HGPA partitioning at K = 2: 2 clusters at ANMI = 0.00329903476904425"
## [1] "HGPA partitioning at K = 3: 3 clusters at ANMI = 0.278691668779803"
## [1] "HGPA partitioning at K = 4: 4 clusters at ANMI = 0.00392992505505839"
## [1] "HGPA partitioning at K = 5: 5 clusters at ANMI = 0.552234460801785"
## [1] "MCLA partitioning at K = 2: 2 clusters at ANMI = 0.568294023177534"
## [1] "MCLA partitioning at K = 3: 3 clusters at ANMI = 0.929094923585274"
## [1] "MCLA partitioning at K = 4: 4 clusters at ANMI = 0.872601957447147"
## [1] "MCLA partitioning at K = 5: 4 clusters at ANMI = 0.923346490477427"
## [1] "CSPA partitioning at K = 2: 2 clusters at ANMI = 0.53144399728197"
## [1] "CSPA partitioning at K = 3: 3 clusters at ANMI = 0.850151780486274"
## [1] "CSPA partitioning at K = 4: 4 clusters at ANMI = 0.665510270422344"
## [1] "CSPA partitioning at K = 5: 5 clusters at ANMI = 0.666022118059772"
## [1] "Optimal number of clusters is 3 with ANMI = 0.929094923585274"

Function SAFE will output a list for Average Normalized Mutual Information (ANMI) metric (Strehl and Ghosh Proceedings of AAAI 2002, Edmonto, Canada, 2002) between each ensemble solution and the individual solutions. The optimal clustering ensemble is selected from the ensemble solution with the highest ANMI value.

cluster.ensemble$Summary
## [1] "Optimal number of clusters is 3 with ANMI = 0.944409180075142"
cluster.ensemble$MCLA[1:10]
##  [1] 1 1 1 1 1 1 1 1 1 1
cluster.ensemble$MCLA_optimal_k
## [1] 3

We can compare the clustering results to the true labels using the Adjusted Rand Index (ARI)

library(cidr)

# Cell labels of ground truth
head(data_SAFE$Zheng.celltype)
## [1] cd56_NK cd56_NK cd56_NK cd56_NK cd56_NK cd56_NK
## Levels: bcell cd56_NK regulatory T
# Calculating ARI for cluster ensemble
adjustedRandIndex(cluster.ensemble$optimal_clustering, data_SAFE$Zheng.celltype)
## [1] 0.9941685

4 Biase dataset

4.1 Setup the input expression matrix

dim(data_SAFE$Biase.expr.expr)
## NULL
data_SAFE$Biase.expr[1:5, 1:5]
##                    GSM1377859 GSM1377860 GSM1377861 GSM1377862 GSM1377863
## ENSMUSG00000000001    25.8078    36.7561    8.87692    24.5712    31.2255
## ENSMUSG00000000028    93.4291    92.1165   94.59080   107.0380   121.4490
## ENSMUSG00000000031     0.0000     0.0000    0.00000     0.0000     0.0000
## ENSMUSG00000000037    37.9544    22.4305   23.34200    42.2728    23.8579
## ENSMUSG00000000049     0.0000     0.0000    0.00000     0.0000     0.0000

4.2 Perform individual clustering

Here we perform single-cell clustering using four popular methods, SC3, CIDR, Seurat and t-SNE + k-means, without filtering any genes or cells. Since there are only 49 cells in Biase dataset, the resolution parameter is set to 1.2 according to our benchmarking results.

cluster.results <- individual_clustering(inputTags = data_SAFE$Biase.expr, datatype = "FPKM",  mt_filter = FALSE, nGene_filter = FALSE, SC3 = TRUE, gene_filter = FALSE, CIDR = TRUE, nPC.cidr = NULL, Seurat = TRUE, nPC.seurat = NULL, seurat_min_cell = 200, resolution_min = 1.2, tSNE = TRUE, dimensions = 3, tsne_min_cells = 200, tsne_min_perplexity = 10, SEED = 123)
## Performing SC3 clustering...
## Estimating k...
## Setting SC3 parameters...
## Calculating distances between the cells...
## Performing transformations and calculating eigenvectors...
## Performing k-means clustering...
## Calculating consensus matrix...
## Performing CIDR clustering...
## Performing Seurat clustering...
## Regressing out: nUMI
## Scaling data matrix
## 1 singletons identified. 3 final clusters.
## Performing tSNE + k-means clustering...

4.3 Cluster ensemble

Using the clustering results, we perform cluster ensemble using all the three partitioning algorithms MCLA, HGPA and CSPA.

cluster.ensemble <- SAFE(cluster_results = cluster.results, program.dir = "~/Documents/single_cell_clustering", MCLA = TRUE, CSPA = TRUE, HGPA = TRUE, SEED = 123)

Here is the list of ANMI results for esemble solution of each K and each partitioning algorithm.

## [1] "HGPA partitioning at K = 2: 2 clusters at ANMI = 0.156896209024547"
## [1] "HGPA partitioning at K = 3: 3 clusters at ANMI = 0.59768416631598"
## [1] "HGPA partitioning at K = 4: 4 clusters at ANMI = 0.614176459706577"
## [1] "MCLA partitioning at K = 2: 2 clusters at ANMI = 0.784102309002763"
## [1] "MCLA partitioning at K = 3: 3 clusters at ANMI = 0.970539568368452"
## [1] "MCLA partitioning at K = 4: 4 clusters at ANMI = 0.971666531448806"
## [1] "CSPA partitioning at K = 2: 2 clusters at ANMI = 0.601004834004939"
## [1] "CSPA partitioning at K = 3: 3 clusters at ANMI = 0.622097187347639"
## [1] "CSPA partitioning at K = 4: 4 clusters at ANMI = 0.590251500201678"
## [1] "Optimal number of clusters is 4 with ANMI = 0.971666531448806"
cluster.ensemble$Summary
## [1] "Optimal number of clusters is 4 with ANMI = 0.971666531448806"
cluster.ensemble$MCLA[1:10]
##  [1] 4 4 4 4 4 4 2 4 4 1
cluster.ensemble$MCLA_optimal_k
## [1] 4

We can compare the clustering results to the true labels using the Adjusted Rand Index (ARI)

library(cidr)

# Cell labels of ground truth
head(data_SAFE$Biase.celltype)
## [1] zygote zygote zygote zygote zygote zygote
## Levels: Four-cell Two-cell zygote
# Calculating ARI for cluster ensemble
adjustedRandIndex(cluster.ensemble$optimal_clustering, data_SAFE$Biase.celltype)
## [1] 0.9850564