Method comparison
Christoph Hafemeister
2020-12-15
Introduction
In sctransform::vst we support several methods to estimate the parameters of the per-gene linear models. Here we briefly go over the methods and compare their results and runtime.
Methods
poisson- does poisson regression and the negative binomial theta parameter is estimated using the response residuals. By default theta is estimated usingMASS::theta.ml, butMASS::theta.mmcan also be used by changing thetheta_estimation_funparameter.qpoisson- does quasi poisson regression to obtain coefficients and overdispersion (phi) and theta is estimated based on phi and the mean fitted value - this is currently the fastest method with results very similar toglmGamPoinb_fast- coefficients and theta are estimated as in thepoissonmethod, but coefficients are then re-estimated using a proper negative binomial model in a call toglmwithfamily = MASS::negative.binomial(theta = theta).nb- coefficients and theta are estimated byMASS::glm.nb.glmGamPoi- coefficients and theta are estimated byglmGamPoi::glm_gp.offset- no regression parameters are learned, but instead an offset model is assumed. The latent variable is set to log_umi and a fixed slope oflog(10)is used (offset). The intercept is given bylog(gene_mean) - log(avg_cell_umi). Theta is set to 100 by default, but can be changed using the theta_given parameteroffset_shared_theta_estimate- like offset above, but the 250 most highly expressed genes with detectio rate of at least 0.5 are used to estimate a theta that is then shared across all genes. Thetas are estimated per individual gene using 5000 randomly selected cells. The final theta used for all genes is then the average.
Comparison
We are going to process a PBMC dataset with all the methods listed above. The dataset is available from 10x Genomics.
cm <- Seurat::Read10X_h5(file = "~/Projects/data_warehouse/raw_public_10x/pbmc_10k_v3_filtered_feature_bc_matrix.h5")
# downsample to speed up compilation of this vignette
set.seed(42)
cm <- cm[, sample(x = ncol(cm), size = 5000)]
message(nrow(cm), " genes across ", ncol(cm), " cells")
#> 33538 genes across 5000 cellsresult_lst <- list()
estimation_methods <- c("poisson", "qpoisson", "nb_fast", "nb", "glmGamPoi", "offset",
"offset_shared_theta_estimate")
for (estimation_method in estimation_methods) {
if (estimation_method %in% c("poisson", "nb_fast")) {
for (theta_estimation_fun in c("theta.ml", "theta.mm")) {
method_name <- paste(estimation_method, theta_estimation_fun, sep = "-")
message(method_name)
set.seed(33)
vst_out <- vst(umi = cm, method = estimation_method, theta_estimation_fun = theta_estimation_fun,
verbosity = 0)
vst_out$y <- NULL
result_lst[[method_name]] <- vst_out
}
} else {
method_name <- estimation_method
message(method_name)
set.seed(33)
vst_out <- vst(umi = cm, method = estimation_method, verbosity = 0)
vst_out$y <- NULL
result_lst[[method_name]] <- vst_out
}
}
#> poisson-theta.ml
#> poisson-theta.mm
#> qpoisson
#> nb_fast-theta.ml
#> nb_fast-theta.mm
#> nb
#> glmGamPoi
#> offset
#> offset_shared_theta_estimateShow how the residual variances compare between methods (plot below shows log10-transformed values)
mat <- sapply(result_lst, function(x) x$gene_attr$residual_variance)
colnames(mat) <- names(result_lst)
ggpairs(data.frame(log10(mat)), progress = FALSE)
Show how the residual means compare between methods
mat <- sapply(result_lst, function(x) x$gene_attr$residual_mean)
colnames(mat) <- names(result_lst)
ggpairs(data.frame(mat), progress = FALSE)
Show the model parameters for all methods
plot_lst <- lapply(names(result_lst), function(method_name) {
plot_model_pars(result_lst[[method_name]], show_theta = TRUE) + ggtitle(method_name)
})
plot(cowplot::plot_grid(plotlist = plot_lst, ncol = 2))
Overall, the regularized parameters and the resulting residuals are very similar across the methods. It is unlikely that the small differences we see above would lead to big differences in downstream analyses like dimensionality reduction and clustering.
How do the methods compare with respect to runtime?
tmp <- lapply(names(result_lst), function(method_name) {
times <- result_lst[[method_name]]$times
delta_t_total <- as.numeric(times$done) - as.numeric(times$start_time)
delta_t_model <- as.numeric(times$reg_model_pars) - as.numeric(times$get_model_pars)
data.frame(method_name, delta_t_total, delta_t_model)
})
df <- do.call(rbind, tmp) %>% mutate(delta_t_rest = delta_t_total - delta_t_model)
df$method_name <- factor(df$method_name, levels = df$method_name, ordered = TRUE)
df <- arrange(df, delta_t_total) %>% mutate(method_name = factor(method_name, unique(method_name)))
melt(df, id.vars = "method_name", measure.vars = c("delta_t_model", "delta_t_rest")) %>%
ggplot(aes(method_name, value, fill = variable)) + geom_bar(position = "stack",
stat = "identity") + theme(axis.text.x = element_text(angle = 45, vjust = 1,
hjust = 1)) + xlab("Method") + ylab("Wall clock time in seconds") + scale_fill_discrete(name = "Part of algorithm",
labels = c("Model fitting", "Rest"))
df
#> method_name delta_t_total delta_t_model delta_t_rest
#> 1 offset 12.77798 0.3969212 12.38106
#> 2 offset_shared_theta_estimate 14.48155 2.1186938 12.36286
#> 3 qpoisson 16.43803 3.7639441 12.67408
#> 4 poisson-theta.mm 19.81587 7.1540921 12.66178
#> 5 glmGamPoi 32.08939 19.0863469 13.00304
#> 6 poisson-theta.ml 70.75487 58.1091411 12.64573
#> 7 nb_fast-theta.mm 72.42107 58.9494340 13.47164
#> 8 nb_fast-theta.ml 124.26923 110.7357450 13.53348
#> 9 nb 642.00292 628.8870690 13.11585The runtime comparison above is not an in-depth benchmark, but it gives a good idea about the relative speed of the different methods. All times above were obtained using a single CPU core. However, the model fitting step supports the Future API for parallel processing.