Contents

1 Introduction

In this vignette we describe the basic usage of the DuoClustering2018 package: how to retrieve data sets and clustering results, and how to construct various plots summarizing the performance of different methods across several data sets.

2 Load the necessary packages

suppressPackageStartupMessages({
  library(ExperimentHub)
  library(SingleCellExperiment)
  library(DuoClustering2018)
  library(plyr)
})

3 Retrieve a data set

The clustering evaluation (Duò, Robinson, and Soneson 2018) is based on 12 data sets (9 real and 3 simulated), which are all provided via ExperimentHub and retrievable via this package. We include the full data sets (after quality filtering of cells and removal of genes with zero counts across all cells) as well as three filtered versions of each data set (by expression, variability and dropout pattern, respectively), each containing 10% of the genes in the full data set.

To get an overview, we can list all records from this package that are available in ExperimentHub:

eh <- ExperimentHub()
## snapshotDate(): 2018-10-31
query(eh, "DuoClustering2018")
## ExperimentHub with 122 records
## # snapshotDate(): 2018-10-31 
## # $dataprovider: Robinson group (UZH), 10x Genomics, Zheng et al (2017), ...
## # $species: Homo sapiens, Mus musculus, NA
## # $rdataclass: data.frame, SingleCellExperiment, list
## # additional mcols(): taxonomyid, genome, description,
## #   coordinate_1_based, maintainer, rdatadateadded, preparerclass,
## #   tags, rdatapath, sourceurl, sourcetype 
## # retrieve records with, e.g., 'object[["EH1499"]]' 
## 
##            title                                               
##   EH1499 | duo_clustering_all_parameter_settings_v1            
##   EH1500 | sce_full_Koh                                        
##   EH1501 | sce_filteredExpr10_Koh                              
##   EH1502 | sce_filteredHVG10_Koh                               
##   EH1503 | sce_filteredM3Drop10_Koh                            
##   ...      ...                                                 
##   EH1651 | clustering_summary_filteredHVG10_SimKumar4hard_v2   
##   EH1652 | clustering_summary_filteredM3Drop10_SimKumar4hard_v2
##   EH1653 | clustering_summary_filteredExpr10_SimKumar8hard_v2  
##   EH1654 | clustering_summary_filteredHVG10_SimKumar8hard_v2   
##   EH1655 | clustering_summary_filteredM3Drop10_SimKumar8hard_v2

The records with names starting in sce_ represent (filtered or unfiltered) data sets (in SingleCellExperiment format). The records with names starting in clustering_summary_ correspond to data.frame objects with clustering results for each of the filtered data sets. Finally, the duo_clustering_all_parameter_settings object contains the parameter settings we used for all the clustering methods. For clustering summaries and parameter settings, the version number (e.g., _v2) corresponds to the version of the publication.

The records can be retrieved using their ExperimentHub ID, e.g.:

eh[["EH1500"]]
## see ?DuoClustering2018 and browseVignettes('DuoClustering2018') for documentation
## downloading 0 resources
## loading from cache 
##     '/home/biocbuild//.ExperimentHub/1500'
## class: SingleCellExperiment 
## dim: 48981 531 
## metadata(1): log.exprs.offset
## assays(3): counts logcounts normcounts
## rownames(48981): ENSG00000000003.14 ENSG00000000005.5 ...
##   ENSG00000283122.1 ENSG00000283124.1
## rowData names(8): is_feature_control mean_counts ... total_counts
##   log10_total_counts
## colnames(531): SRR3952323 SRR3952325 ... SRR3952970 SRR3952971
## colData names(14): Run LibraryName ... libsize.drop feature.drop
## reducedDimNames(2): PCA TSNE
## spikeNames(0):

Alternatively, the shortcut functions provided by this package can be used:

sce_filteredExpr10_Koh()
## snapshotDate(): 2018-10-31
## see ?DuoClustering2018 and browseVignettes('DuoClustering2018') for documentation
## downloading 0 resources
## loading from cache 
##     '/home/biocbuild//.ExperimentHub/1501'
## class: SingleCellExperiment 
## dim: 4898 531 
## metadata(1): log.exprs.offset
## assays(3): counts logcounts normcounts
## rownames(4898): ENSG00000198804.2 ENSG00000210082.2 ...
##   ENSG00000072134.15 ENSG00000090061.17
## rowData names(8): is_feature_control mean_counts ... total_counts
##   log10_total_counts
## colnames(531): SRR3952323 SRR3952325 ... SRR3952970 SRR3952971
## colData names(14): Run LibraryName ... pct_counts_top_500_features
##   is_cell_control
## reducedDimNames(2): PCA TSNE
## spikeNames(0):

4 Read a set of clustering results

For each included data set, we have applied a range of clustering methods (see the run_clustering vignette for more details on how this was done, and how to apply additional methods). As mentioned above, the results of these clusterings are also available from ExperimentHub, and can be loaded either by their ExperimentHub ID or using the provided shortcut functions, as above. For simplicity, the results of all methods for a given data set are combined into a single object. As an illustration, we load the clustering summaries for two different data sets (Koh and Zhengmix4eq), each with two different gene filterings (Expr10 and HVG10):

res <- plyr::rbind.fill(
  clustering_summary_filteredExpr10_Koh_v2(),
  clustering_summary_filteredHVG10_Koh_v2(),
  clustering_summary_filteredExpr10_Zhengmix4eq_v2(),
  clustering_summary_filteredHVG10_Zhengmix4eq_v2()
)
## snapshotDate(): 2018-10-31
## see ?DuoClustering2018 and browseVignettes('DuoClustering2018') for documentation
## downloading 0 resources
## loading from cache 
##     '/home/biocbuild//.ExperimentHub/1620'
## snapshotDate(): 2018-10-31
## see ?DuoClustering2018 and browseVignettes('DuoClustering2018') for documentation
## downloading 0 resources
## loading from cache 
##     '/home/biocbuild//.ExperimentHub/1621'
## snapshotDate(): 2018-10-31
## see ?DuoClustering2018 and browseVignettes('DuoClustering2018') for documentation
## downloading 0 resources
## loading from cache 
##     '/home/biocbuild//.ExperimentHub/1638'
## snapshotDate(): 2018-10-31
## see ?DuoClustering2018 and browseVignettes('DuoClustering2018') for documentation
## downloading 0 resources
## loading from cache 
##     '/home/biocbuild//.ExperimentHub/1639'
dim(res)
## [1] 5625885      10

The resulting data.frame contains 10 columns:

head(res)
##                  dataset    method       cell run k resolution cluster
## 1 sce_filteredExpr10_Koh PCAKmeans SRR3952323   1 2         NA       1
## 2 sce_filteredExpr10_Koh PCAKmeans SRR3952325   1 2         NA       1
## 3 sce_filteredExpr10_Koh PCAKmeans SRR3952326   1 2         NA       1
## 4 sce_filteredExpr10_Koh PCAKmeans SRR3952327   1 2         NA       1
## 5 sce_filteredExpr10_Koh PCAKmeans SRR3952328   1 2         NA       1
## 6 sce_filteredExpr10_Koh PCAKmeans SRR3952329   1 2         NA       1
##   trueclass est_k elapsed
## 1    H7hESC    NA  14.318
## 2    H7hESC    NA  14.318
## 3    H7hESC    NA  14.318
## 4    H7hESC    NA  14.318
## 5    H7hESC    NA  14.318
## 6    H7hESC    NA  14.318

5 Define consistent method colors

For some of the plots generated below, the points will be colored according to the clustering method. We can enforce a consistent set of colors for the methods by defining a named vector of colors to use for all plots.

method_colors <- c(CIDR = "#332288", FlowSOM = "#6699CC", PCAHC = "#88CCEE", 
            PCAKmeans = "#44AA99", pcaReduce = "#117733",
            RtsneKmeans = "#999933", Seurat = "#DDCC77", SC3svm = "#661100", 
            SC3 = "#CC6677", TSCAN = "grey34", ascend = "orange", SAFE = "black",
            monocle = "red", RaceID2 = "blue")

6 Plot

Each plotting function described below returns a list of ggplot objects. These can be plotted directly, or further modified if desired.

6.1 Performance

The plot_performance() function generates plots related to the performance of the clustering methods. We quantify performance using the adjusted Rand Index (ARI) (Hubert and Arabie 1985), comparing the obtained clustering to the true clusters. As we noted in the publication (Duò, Robinson, and Soneson 2018), defining a true partitioning of the cells is difficult, since they can often be grouped together in several different, but still interpretable, ways. We refer to our paper for more information on how the true clusters were defined for each of the data sets.

perf <- plot_performance(res, method_colors = method_colors)
names(perf)
## [1] "median_ari_vs_k"           "scatter_time_vs_ari_truek"
## [3] "median_ari_heatmap_truek"  "median_ari_heatmap_bestk" 
## [5] "median_ari_heatmap_estk"
perf$median_ari_vs_k
## Warning: Removed 4 rows containing missing values (geom_path).

perf$median_ari_heatmap_truek

6.2 Stability

The plot_stability() function evaluates the stability of the clustering results from each method, with respect to random starts. Each method was run five times on each data set (for each k), and we quantify the stability by comparing each pair of such runs using the adjusted Rand Index.

stab <- plot_stability(res, method_colors = method_colors)
## Aggregation function missing: defaulting to length
## Aggregation function missing: defaulting to length
## Aggregation function missing: defaulting to length
## Aggregation function missing: defaulting to length
## Aggregation function missing: defaulting to length
## Aggregation function missing: defaulting to length
## Aggregation function missing: defaulting to length
## Aggregation function missing: defaulting to length
## Aggregation function missing: defaulting to length
## Aggregation function missing: defaulting to length
## Aggregation function missing: defaulting to length
## Aggregation function missing: defaulting to length
## Aggregation function missing: defaulting to length
## Aggregation function missing: defaulting to length
names(stab)
## [1] "stability_vs_k"          "stability_truek"        
## [3] "stability_heatmap_truek"
stab$stability_vs_k
## Warning: Removed 4 rows containing missing values (geom_path).

stab$stability_heatmap_truek

6.3 Entropy

In order to evaluate the tendency of the clustering methods to favor equally sized clusters, we calculate the Shannon entropy (Shannon 1948) of each clustering solution (based on the proportions of cells in the different clusters) and plot this using the plot_entropy() function. Since the maximal entropy that can be obtained depends on the number of clusters, we use normalized entropies, defined by dividing the observed entropy by log2(k). We also compare the entropies for the clusterings to the entropy of the true partition for each data set.

entr <- plot_entropy(res, method_colors = method_colors)
names(entr)
## [1] "entropy_vs_k"           "entropy_vs_ari"        
## [3] "normentropy"            "deltaentropy_truek"    
## [5] "deltanormentropy_truek"
entr$entropy_vs_k

entr$normentropy

6.4 Timing

The plot_timing() function plots various aspects of the timing of the different methods.

timing <- plot_timing(res, method_colors = method_colors, 
                      scaleMethod = "RtsneKmeans")
names(timing)
## [1] "time_boxplot"           "time_normalized_by_ref"
## [3] "time_vs_k"
timing$time_normalized_by_ref
## Warning: Removed 32 rows containing non-finite values (stat_boxplot).

6.5 Differences in k

Most performance evaluations above are performed on the clustering solutions obtained by imposing the “true” number of clusters. The plot_k_diff() function evaluates the difference between the true number of cluster and the number of clusters giving the best agreement with the true partition, as well as the difference between the estimated and the true number of clusters, for the methods that allow estimation of k.

kdiff <- plot_k_diff(res, method_colors = method_colors)
names(kdiff)
## [1] "diff_kmax_ktrue" "diff_kest_ktrue"
kdiff$diff_kest_ktrue
## `stat_bindot()` using `bins = 30`. Pick better value with `binwidth`.

7 Session info

sessionInfo()
## R version 3.5.1 Patched (2018-07-12 r74967)
## Platform: x86_64-pc-linux-gnu (64-bit)
## Running under: Ubuntu 16.04.5 LTS
## 
## Matrix products: default
## BLAS: /home/biocbuild/bbs-3.8-bioc/R/lib/libRblas.so
## LAPACK: /home/biocbuild/bbs-3.8-bioc/R/lib/libRlapack.so
## 
## locale:
##  [1] LC_CTYPE=en_US.UTF-8       LC_NUMERIC=C              
##  [3] LC_TIME=en_US.UTF-8        LC_COLLATE=C              
##  [5] LC_MONETARY=en_US.UTF-8    LC_MESSAGES=en_US.UTF-8   
##  [7] LC_PAPER=en_US.UTF-8       LC_NAME=C                 
##  [9] LC_ADDRESS=C               LC_TELEPHONE=C            
## [11] LC_MEASUREMENT=en_US.UTF-8 LC_IDENTIFICATION=C       
## 
## attached base packages:
## [1] parallel  stats4    stats     graphics  grDevices utils     datasets 
## [8] methods   base     
## 
## other attached packages:
##  [1] plyr_1.8.4                  ExperimentHub_1.8.0        
##  [3] AnnotationHub_2.14.0        bindrcpp_0.2.2             
##  [5] tidyr_0.8.2                 dplyr_0.7.7                
##  [7] DuoClustering2018_1.0.0     SingleCellExperiment_1.4.0 
##  [9] SummarizedExperiment_1.12.0 DelayedArray_0.8.0         
## [11] BiocParallel_1.16.0         matrixStats_0.54.0         
## [13] Biobase_2.42.0              GenomicRanges_1.34.0       
## [15] GenomeInfoDb_1.18.0         IRanges_2.16.0             
## [17] S4Vectors_0.20.0            BiocGenerics_0.28.0        
## [19] BiocStyle_2.10.0           
## 
## loaded via a namespace (and not attached):
##  [1] viridis_0.5.1                 httr_1.3.1                   
##  [3] bit64_0.9-7                   viridisLite_0.3.0            
##  [5] shiny_1.1.0                   assertthat_0.2.0             
##  [7] interactiveDisplayBase_1.20.0 BiocManager_1.30.3           
##  [9] blob_1.1.1                    GenomeInfoDbData_1.2.0       
## [11] yaml_2.2.0                    pillar_1.3.0                 
## [13] RSQLite_2.1.1                 backports_1.1.2              
## [15] lattice_0.20-35               glue_1.3.0                   
## [17] digest_0.6.18                 promises_1.0.1               
## [19] XVector_0.22.0                colorspace_1.3-2             
## [21] htmltools_0.3.6               httpuv_1.4.5                 
## [23] Matrix_1.2-14                 pkgconfig_2.0.2              
## [25] bookdown_0.7                  zlibbioc_1.28.0              
## [27] xtable_1.8-3                  purrr_0.2.5                  
## [29] scales_1.0.0                  later_0.7.5                  
## [31] tibble_1.4.2                  ggplot2_3.1.0                
## [33] lazyeval_0.2.1                mime_0.6                     
## [35] magrittr_1.5                  crayon_1.3.4                 
## [37] mclust_5.4.1                  memoise_1.1.0                
## [39] evaluate_0.12                 ggthemes_4.0.1               
## [41] tools_3.5.1                   stringr_1.3.1                
## [43] munsell_0.5.0                 AnnotationDbi_1.44.0         
## [45] compiler_3.5.1                rlang_0.3.0.1                
## [47] grid_3.5.1                    RCurl_1.95-4.11              
## [49] labeling_0.3                  bitops_1.0-6                 
## [51] rmarkdown_1.10                gtable_0.2.0                 
## [53] curl_3.2                      DBI_1.0.0                    
## [55] reshape2_1.4.3                R6_2.3.0                     
## [57] gridExtra_2.3                 knitr_1.20                   
## [59] bit_1.1-14                    bindr_0.1.1                  
## [61] rprojroot_1.3-2               stringi_1.2.4                
## [63] Rcpp_0.12.19                  tidyselect_0.2.5             
## [65] xfun_0.4

References

Duò, A, MD Robinson, and D Soneson. 2018. “A Systematic Performance Evaluation of Clustering Methods for Single-Cell RNA-seq Data.” F1000Research 7:1141.

Hubert, L, and P Arabie. 1985. “Comparing Partitions.” Journal of Classification 2 (1):193–218.

Shannon, C E. 1948. “A Mathematical Theory of Communication.” The Bell System Technical Journal 27 (3):379–423.