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(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
camp_dat <- dataset
#Remove the "unknown" cell type
to_keep <- which(!(colData(camp_dat)$cell_type1 %in% c("Unknown")))
camp_dat <- camp_dat[, to_keep]
#Creates FPKM assay
logs <- (as.matrix(assays(camp_dat)[["logcounts"]]))
fpkm(camp_dat) <- apply(logs, c(1,2), function(x) 2^x-1)
#Find the genes with all zero entries and remove
keep_features <- rowSums(fpkm(camp_dat) > 0) > 0
camp_dat <- camp_dat[keep_features, ]
#Use Scran normalisation
computeSumFactors(camp_dat, assay.type = "fpkm")
## class: SingleCellExperiment
## dim: 16450 553
## metadata(0):
## assays(2): logcounts fpkm
## rownames(16450): A1BG A1CF ... ZZEF1 ZZZ3
## rowData names(1): feature_symbol
## colnames(553): 12wpc_A1_c1_fetal 12wpc_A1_c2_fetal ... H8_hOrg37d
## H8_hOrg41d
## colData names(5): Species cell_type1 Source age batch
## reducedDimNames(0):
## spikeNames(1): ERCC
normalize(camp_dat, exprs_values = "fpkm")
## Warning in .local(object, ...): using library sizes as size factors
## Warning in .get_all_sf_sets(object): spike-in set 'ERCC' should have its
## own size factors
## class: SingleCellExperiment
## dim: 16450 553
## metadata(1): log.exprs.offset
## assays(2): logcounts fpkm
## rownames(16450): A1BG A1CF ... ZZEF1 ZZZ3
## rowData names(1): feature_symbol
## colnames(553): 12wpc_A1_c1_fetal 12wpc_A1_c2_fetal ... H8_hOrg37d
## H8_hOrg41d
## colData names(5): Species cell_type1 Source age batch
## reducedDimNames(0):
## spikeNames(1): ERCC
#Remove the spikes
is.spike <-grepl("^ERCC", rownames(camp_dat))
table(is.spike)
## is.spike
## FALSE TRUE
## 16440 10
camp_dat <- camp_dat[!is.spike, ]
#Start the time here
start_time <- Sys.time()
#Extract the needed variables for scGPS
camp_dat_exprs <- assays(camp_dat)[["logcounts"]]
camp_dat_cellnames <- colnames(camp_dat)
camp_dat_cellnames <- data.frame("cellBarcodes" = camp_dat_cellnames)
camp_dat_GeneMetaData <- rownames(camp_dat)
camp_dat_GeneMetaData <- data.frame("GeneSymbol" = camp_dat_GeneMetaData)
#Store Data in scGPS format
mixedpop <- new_summarized_scGPS_object(ExpressionMatrix = camp_dat_exprs, GeneMetadata = camp_dat_GeneMetaData, CellMetadata = camp_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
## 553 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
## 57178NANA
## 3
## 57178278NA
## 4
## 57178278434
## 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
#Make a dataframe with the results we want to examine
cell_types1 <- colData(camp_dat)$cell_type1
label_list <- unlist(camp_dat_cellnames$cellBarcodes)
cluster_list <- unlist(CORE_cluster_bagging$Cluster[CORE_cluster_bagging$optimal_index])
compare_frame <- data.frame("Gene_label" = label_list, "type1" = cell_types1,"cluster" = cluster_list)
#Find the Adjusted Rand Index
AdjustedRandIndex_SCORE <- mclust::adjustedRandIndex(compare_frame$type1, compare_frame$cluster)
HighResRand <- mclust::adjustedRandIndex(compare_frame$type1, unlist(CORE_cluster_bagging$Cluster[1]))
#Store the estimated k from the bagging runs
estimated_k_SCORE <- CORE_cluster_bagging$optimalMax
#Remove unwanted data
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.5589109
## --- dataset ---
## class: SingleCellExperiment
## dim: 18927 734
## metadata(0):
## assays(1): logcounts
## rownames(18927): A1BG A1CF ... ZZEF1 ZZZ3
## rowData names(1): feature_symbol
## colnames(734): 12wpc_A1_c1_fetal 12wpc_A1_c2_fetal ... H8_hOrg41d
## H9_hOrg65d
## colData names(5): Species cell_type1 Source age batch
## reducedDimNames(0):
## spikeNames(1): ERCC
## --- estimated_k_SCORE ---
## [1] 4
## --- HighResRand ---
## [1] 0.5972624
## --- time_difference_SCORE ---
## Time difference of 32.61574 secs
SC3
#Load everything for SC3
library(SC3)
library(scater)
#Retrieve the dataset
sce <- dataset
#Remove the "unknown" cell type
to_keep <- which(!(colData(sce)$cell_type1 %in% c("Unknown")))
sce <- sce[, to_keep]
#Creates FPKM assay
logs <- (as.matrix(assays(sce)[["logcounts"]]))
fpkm(sce) <- apply(logs, c(1,2), function(x) 2^x-1)
counts(sce) <- fpkm(sce)
#Find the genes with all zero entries and remove
keep_features <- rowSums(fpkm(sce) > 0) > 0
sce <- sce[keep_features, ]
#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
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
cell_types1 <- colData(sce)$cell_type1
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, "type1" = cell_types1, "cluster" = cluster_list)
#Find the Adjusted Rand Index
AdjustedRandIndex_SC3 <- mclust::adjustedRandIndex(compare_frame$type1, 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.5556678
## --- dataset ---
## class: SingleCellExperiment
## dim: 18927 734
## metadata(0):
## assays(1): logcounts
## rownames(18927): A1BG A1CF ... ZZEF1 ZZZ3
## rowData names(1): feature_symbol
## colnames(734): 12wpc_A1_c1_fetal 12wpc_A1_c2_fetal ... H8_hOrg41d
## H9_hOrg65d
## colData names(5): Species cell_type1 Source age batch
## reducedDimNames(0):
## spikeNames(1): ERCC
## --- SC3_k_estimate ---
## [1] 10
## --- time_difference_SC3 ---
## Time difference of 2.693245 mins