Benchmaking for the koh dataset

Load the dataset to use

dataset <- readRDS(url("http://imlspenticton.uzh.ch/robinson_lab/conquer/data-mae/SRP073808.rds"))

scGPS

#Load everyting for scGPS Benchmarking
library(scGPS)
## Loading required package: SummarizedExperiment
## Loading required package: GenomicRanges
## Loading required package: stats4
## Loading required package: BiocGenerics
## Loading required package: parallel
## 
## Attaching package: 'BiocGenerics'
## The following objects are masked from 'package:parallel':
## 
##     clusterApply, clusterApplyLB, clusterCall, clusterEvalQ,
##     clusterExport, clusterMap, parApply, parCapply, parLapply,
##     parLapplyLB, parRapply, parSapply, parSapplyLB
## The following objects are masked from 'package:stats':
## 
##     IQR, mad, sd, var, xtabs
## The following objects are masked from 'package:base':
## 
##     anyDuplicated, append, as.data.frame, basename, cbind,
##     colnames, dirname, do.call, duplicated, eval, evalq, Filter,
##     Find, get, grep, grepl, intersect, is.unsorted, lapply, Map,
##     mapply, match, mget, order, paste, pmax, pmax.int, pmin,
##     pmin.int, Position, rank, rbind, Reduce, rownames, sapply,
##     setdiff, sort, table, tapply, union, unique, unsplit, which,
##     which.max, which.min
## Loading required package: S4Vectors
## 
## Attaching package: 'S4Vectors'
## The following object is masked from 'package:base':
## 
##     expand.grid
## Loading required package: IRanges
## 
## Attaching package: 'IRanges'
## The following object is masked from 'package:grDevices':
## 
##     windows
## Loading required package: GenomeInfoDb
## Loading required package: Biobase
## Welcome to Bioconductor
## 
##     Vignettes contain introductory material; view with
##     'browseVignettes()'. To cite Bioconductor, see
##     'citation("Biobase")', and for packages 'citation("pkgname")'.
## Loading required package: DelayedArray
## Loading required package: matrixStats
## 
## Attaching package: 'matrixStats'
## The following objects are masked from 'package:Biobase':
## 
##     anyMissing, rowMedians
## Loading required package: BiocParallel
## 
## Attaching package: 'DelayedArray'
## The following objects are masked from 'package:matrixStats':
## 
##     colMaxs, colMins, colRanges, rowMaxs, rowMins, rowRanges
## The following objects are masked from 'package:base':
## 
##     aperm, apply, rowsum
## Loading required package: dynamicTreeCut
## Loading required package: SingleCellExperiment
## Registered S3 methods overwritten by 'ggplot2':
##   method         from 
##   [.quosures     rlang
##   c.quosures     rlang
##   print.quosures rlang
library("MultiAssayExperiment")
library("scater")
## Loading required package: ggplot2
## 
## Attaching package: 'scater'
## The following object is masked from 'package:S4Vectors':
## 
##     rename
## The following object is masked from 'package:stats':
## 
##     filter
library("scran")

#Retrieve the dataset
koh_dat <- dataset

#Exract the gene-level length-scaled TPMs
koh_expr <- assays(experiments(koh_dat)[["gene"]])[["count_lstpm"]]

#Extract the phenotype data.
phn <- colData(koh_dat)
phn$phenoid <- as.character(phn$LibraryName)
table(phn$phenoid)
## 
##           H7_derived_APS         H7_derived_D2LtM  H7_derived_D3GARPpCrdcM 
##                       64                       55                       52 
## H7_derived_D5CntrlDrmmtm      H7_derived_DLL1pPXM          H7_derived_ESMT 
##                       67                       75                       45 
##           H7_derived_MPS        H7_derived_Sclrtm  H7_dreived_D2.25_Smtmrs 
##                       60                       69                       87 
##                   H7hESC 
##                       77
#Create single cell experiment
stopifnot(all(colnames(koh_expr) == rownames(phn)))
SCE <- SingleCellExperiment(
  assays = list(counts = koh_expr),
  colData = phn
)

#Remove features with no gene expression
keep_features <- rowSums(counts(SCE) > 0) > 0
SCE <- SCE[keep_features, ]

#Use scran normalisation
SCE <- computeSumFactors(SCE)
SCE <- normalize(SCE, exprs_values = "counts", return_log = TRUE) 

#Create a count per million assay
cpm(SCE) <- calculateCPM(SCE)

#Remove spikes
is.spike <-grepl("^ERCC", rownames(SCE))
SCE <- SCE[!is.spike, ]

#Start the time here
start_time <- Sys.time()

#Extract the needed variables
koh_dat_exprs <- assays(SCE)[["logcounts"]]
koh_dat_cellnames <- colnames(SCE)
koh_dat_cellnames <- data.frame("cellBarcodes" = koh_dat_cellnames)
koh_dat_GeneMetaData <- rownames(SCE)
koh_dat_GeneMetaData <- data.frame("GeneSymbol" = koh_dat_GeneMetaData)

#Store Data in scGPS format
mixedpop <- new_summarized_scGPS_object(ExpressionMatrix = koh_dat_exprs, GeneMetadata = koh_dat_GeneMetaData, CellMetadata = koh_dat_cellnames)

#Cluster and plot data using SCORE
CORE_cluster_bagging <- CORE_bagging(mixedpop, remove_outlier = c(0), PCA=FALSE)
## Performing 1 round of filtering
## Identifying top variable genes
## Calculating distance matrix
## Performing hierarchical clustering
## Finding clustering information
## No more outliers detected in filtering round 1
## Identifying top variable genes
## Calculating distance matrix
## Performing hierarchical clustering
## Finding clustering information
## 651 cells left after filtering
## Running 20 bagging runs, with 0.8 subsampling...
## Done clustering, moving to stability calculation...
## Done finding optimal clustering
plot_CORE(CORE_cluster_bagging$tree, list_clusters = CORE_cluster_bagging$Cluster)

plot_optimal_CORE(original_tree= CORE_cluster_bagging$tree, optimal_cluster = unlist(CORE_cluster_bagging$Cluster[CORE_cluster_bagging$optimal_index]), shift = -100)
## Ordering and assigning labels...
## 2
## 74173NANANANANA
## 3
## 74173247NANANANA
## 4
## 74173247325NANANA
## 5
## 74173247325399NANA
## 6
## 74173247325399466NA
## 7
## 74173247325399466569
## Plotting the colored dendrogram now....
## Plotting the bar underneath now....

#Stop the time here
end_time <- Sys.time()
time_difference_SCORE <- end_time - start_time

#Find data needed for comparisons and store in data frame
phenoid_list <- unlist(colData(SCE)$phenoid)
label_list <- unlist(koh_dat_cellnames$cellBarcodes)
cluster_list <- unlist(CORE_cluster_bagging$Cluster[CORE_cluster_bagging$optimal_index])
compare_frame <- data.frame("Gene_label" = label_list, "phenoid_list" = phenoid_list, "cluster" = cluster_list)

#Find the adjusted rand index
AdjustedRandIndex_SCORE <- mclust::adjustedRandIndex(compare_frame$phenoid_list, compare_frame$cluster)
HighResRand <- mclust::adjustedRandIndex(compare_frame$phenoid_list, unlist(CORE_cluster_bagging$Cluster[1]))
estimated_k_SCORE <- CORE_cluster_bagging$optimalMax

#Remove clutter from the environment
rm(list = setdiff(ls(), c("AdjustedRandIndex_SCORE", "time_difference_SCORE", "estimated_k_SCORE", "HighResRand", "dataset")))
for ( obj in ls() ) { cat('---',obj,'---\n'); print(get(obj)) }
## --- AdjustedRandIndex_SCORE ---
## [1] 0.5650893
## --- dataset ---
## A MultiAssayExperiment object of 2 listed
##  experiments with user-defined names and respective classes. 
##  Containing an ExperimentList class object of length 2: 
##  [1] gene: RangedSummarizedExperiment with 65218 rows and 651 columns 
##  [2] tx: RangedSummarizedExperiment with 216423 rows and 651 columns 
## Features: 
##  experiments() - obtain the ExperimentList instance 
##  colData() - the primary/phenotype DataFrame 
##  sampleMap() - the sample availability DataFrame 
##  `$`, `[`, `[[` - extract colData columns, subset, or experiment 
##  *Format() - convert into a long or wide DataFrame 
##  assays() - convert ExperimentList to a SimpleList of matrices
## --- estimated_k_SCORE ---
## [1] 7
## --- HighResRand ---
## [1] 0.6612571
## --- time_difference_SCORE ---
## Time difference of 41.76453 secs

SC3

#Load everything for SC3
library("SC3")
library("MultiAssayExperiment")
library("scater")

#Retrieve the dataset
koh_dat <- dataset

#Exract the gene-level length-scaled TPMs
koh_expr <- assays(experiments(koh_dat)[["gene"]])[["count_lstpm"]]

#Extract th phenotype data.
phn <- colData(koh_dat)
phn$phenoid <- as.character(phn$LibraryName)

#Create single cell experiment
stopifnot(all(colnames(koh_expr) == rownames(phn)))
SCE <- SingleCellExperiment(
  assays = list(counts = koh_expr),
  colData = phn
)

#Find the genes with all zero entries and remove
keep_features <- rowSums(counts(SCE) > 0) > 0
SCE <- SCE[keep_features, ]

#Create "logcounts" assays
SCE <- normalize(SCE, exprs_values = "counts", return_log = TRUE)
## Warning in .local(object, ...): using library sizes as size factors
#Remove the spikes
is.spike <-grepl("^ERCC", rownames(SCE))
SCE <- SCE[!is.spike, ]

#Start the time here
start_time <- Sys.time()

#Run sc3 with an estimation for k
rowData(SCE)$feature_symbol <- rownames(counts(SCE))
SCE <- sc3_prepare(SCE, n_cores = 1, gene_filter = TRUE)
## Setting SC3 parameters...
SCE <- sc3_estimate_k(SCE)
## Estimating k...
SC3_k_estimate <- as.integer(unlist(metadata(SCE)$sc3$k_estimation))
SCE <- sc3_calc_dists(SCE)
## Calculating distances between the cells...
SCE <- sc3_calc_transfs(SCE)
## Performing transformations and calculating eigenvectors...
SCE <- sc3_kmeans(SCE, ks = SC3_k_estimate)
## Performing k-means clustering...
SCE <- sc3_calc_consens(SCE)
## Calculating consensus matrix...
#Stop the time here
end_time <- Sys.time()
time_difference_SC3 <- end_time - start_time

#Make a dataframe with the results we want to examine
phenoid_list <- colData(SCE)$phenoid
label_list <- rownames(colData(SCE))
cluster_list <- as.numeric(colData(SCE)[, paste0("sc3_", SC3_k_estimate, "_clusters")])
compare_frame <- data.frame("Gene_label" = label_list, "phenoid_list" = phenoid_list, "cluster" = cluster_list)

#Find the Adjusted Rand Index
AdjustedRandIndex_SC3 <- mclust::adjustedRandIndex(compare_frame$phenoid_list, compare_frame$cluster)

#Remove unwanted data
rm(list = setdiff(ls(), c("AdjustedRandIndex_SC3", "time_difference_SC3", "SC3_k_estimate", "dataset")))
for ( obj in ls() ) { cat('---',obj,'---\n'); print(get(obj)) }
## --- AdjustedRandIndex_SC3 ---
## [1] 0.8239352
## --- dataset ---
## A MultiAssayExperiment object of 2 listed
##  experiments with user-defined names and respective classes. 
##  Containing an ExperimentList class object of length 2: 
##  [1] gene: RangedSummarizedExperiment with 65218 rows and 651 columns 
##  [2] tx: RangedSummarizedExperiment with 216423 rows and 651 columns 
## Features: 
##  experiments() - obtain the ExperimentList instance 
##  colData() - the primary/phenotype DataFrame 
##  sampleMap() - the sample availability DataFrame 
##  `$`, `[`, `[[` - extract colData columns, subset, or experiment 
##  *Format() - convert into a long or wide DataFrame 
##  assays() - convert ExperimentList to a SimpleList of matrices
## --- SC3_k_estimate ---
## [1] 18
## --- time_difference_SC3 ---
## Time difference of 3.238979 mins