UniFrac
Here we’ll compare the time and memory taken by the unweighted,
weighted, weight normalized, generalized, and variance adjusted UniFrac
functions from the abdiv,ampvis2
ecodive, GUniFrac, phyloseq,
phyloregion and picante R packages. Each of
the functions will be run 10 times to time them for speed, and 1 time to
analyze memory usage.
Click to reveal R code.
## Unweighted UniFrac
u_unifrac_res <- rbind(
local({
# cluster for phyloseq
cl <- parallel::makeCluster(n_cpus)
doParallel::registerDoParallel(cl)
on.exit(parallel::stopCluster(cl))
bench::mark(
iterations = 10,
'abdiv' = cleanup(pairwise(abdiv::unweighted_unifrac, hmp50_mtx, hmp50_tree)),
'ecodive' = cleanup(ecodive::unweighted_unifrac(hmp50_mtx, hmp50_tree)),
'GUniFrac' = cleanup(as.dist(GUniFrac::GUniFrac(hmp50_tmtx, hmp50_tree, alpha=1, verbose=FALSE)[[1]][,,2])),
'phyloregion' = cleanup(phyloregion::unifrac(hmp50_tmtx, hmp50_tree)),
'phyloseq' = cleanup(phyloseq::UniFrac(hmp50_phy, weighted=FALSE, normalized=FALSE, parallel=TRUE)),
'picante' = cleanup(picante::unifrac(hmp50_tmtx, hmp50_tree)) )
}),
# ampvis2 conflicts with phyloseq cluster, so run separately
bench::mark(
iterations = 10,
'ampvis2' = {
cleanup(ampvis2:::dist.unifrac(hmp50_mtx, hmp50_tree, weighted=FALSE, normalise=FALSE, num_threads=n_cpus))
doParallel::stopImplicitCluster() } )
)
u_unifrac_res[,1:9]
#> # A tibble: 5 × 13
#> expression min median `itr/sec` mem_alloc `gc/sec` n_itr n_gc total_time
#> <bch:expr> <bch:tm> <bch:tm> <dbl> <bch:byt> <dbl> <int> <dbl> <bch:tm>
#> 1 abdiv 13.84s 14.38s 0.0676 20.1GB 1.87 10 277 2.46m
#> 2 ecodive 5.19ms 5.31ms 184. 770.5KB 0 10 0 54.23ms
#> 3 GUniFrac 77.89ms 80.4ms 11.7 92.1MB 1.17 10 1 858.32ms
#> 4 phyloseq 292.07ms 327.01ms 2.49 49.9MB 0 10 0 4.02s
#> 5 ampvis2 3.36s 3.44s 0.288 49.8MB 0.0320 9 1 31.29s
ggplot(u_unifrac_res, aes(x = median, y = mem_alloc)) +
geom_point() +
geom_label_repel(aes(label = as.character(expression))) +
labs(
title = 'Unweighted UniFrac Implementations',
subtitle = '50 sample all-vs-all benchmarking on six CPU cores',
x = 'Median Calculation Time (log scale; n=10)',
y = 'Memory Allocated\n(log scale)' ) +
theme_bw()
## Weighted UniFrac
w_unifrac_res <- rbind(
local({
# cluster for phyloseq
cl <- parallel::makeCluster(n_cpus)
doParallel::registerDoParallel(cl)
on.exit(parallel::stopCluster(cl))
bench::mark(
iterations = 10,
'abdiv' = cleanup(pairwise(abdiv::weighted_unifrac, hmp50_mtx, hmp50_tree)),
'ecodive' = cleanup(ecodive::weighted_unifrac(hmp50_mtx, hmp50_tree)),
'phyloseq' = cleanup(phyloseq::UniFrac(hmp50_phy, weighted=TRUE, normalized=FALSE, parallel=TRUE)) )
}),
# ampvis2 conflicts with phyloseq cluster, so run separately
bench::mark(
iterations = 10,
'ampvis2' = {
cleanup(ampvis2:::dist.unifrac(hmp50_mtx, hmp50_tree, weighted=TRUE, normalise=FALSE, num_threads=n_cpus))
doParallel::stopImplicitCluster() } )
)
## Weighted Normalized UniFrac
wn_unifrac_res <- rbind(
local({
# cluster for phyloseq
cl <- parallel::makeCluster(n_cpus)
doParallel::registerDoParallel(cl)
on.exit(parallel::stopCluster(cl))
bench::mark(
iterations = 10,
'abdiv' = cleanup(pairwise(abdiv::weighted_normalized_unifrac, hmp50_mtx, hmp50_tree)),
'ecodive' = cleanup(ecodive::weighted_normalized_unifrac(hmp50_mtx, hmp50_tree)),
'GUniFrac' = cleanup(as.dist(GUniFrac::GUniFrac(hmp50_tmtx, hmp50_tree, alpha=1, verbose=FALSE)[[1]][,,1])),
'phyloseq' = cleanup(phyloseq::UniFrac(hmp50_phy, weighted=TRUE, normalized=TRUE, parallel=TRUE)) )
}),
# ampvis2 conflicts with phyloseq cluster, so run separately
bench::mark(
iterations = 10,
'ampvis2' = {
cleanup(ampvis2:::dist.unifrac(hmp50_mtx, hmp50_tree, weighted=TRUE, normalise=TRUE, num_threads=n_cpus))
doParallel::stopImplicitCluster() } )
)
## Weighted Normalized UniFrac
g_unifrac_res <- rbind(
local({
# cluster for phyloseq
cl <- parallel::makeCluster(n_cpus)
doParallel::registerDoParallel(cl)
on.exit(parallel::stopCluster(cl))
bench::mark(
iterations = 10,
'abdiv' = cleanup(pairwise(abdiv::generalized_unifrac, hmp50_mtx, hmp50_tree, alpha=0.5)),
'ecodive' = cleanup(ecodive::generalized_unifrac(hmp50_mtx, hmp50_tree, alpha=0.5)),
'GUniFrac' = cleanup(as.dist(GUniFrac::GUniFrac(hmp50_tmtx, hmp50_tree, alpha=0.5, verbose=FALSE)[[1]][,,1])) )
})
)
## Variance Adjusted UniFrac
va_unifrac_res <- rbind(
local({
# cluster for phyloseq
cl <- parallel::makeCluster(n_cpus)
doParallel::registerDoParallel(cl)
on.exit(parallel::stopCluster(cl))
bench::mark(
iterations = 10,
'abdiv' = cleanup(pairwise(abdiv::variance_adjusted_unifrac, hmp50_mtx, hmp50_tree)),
'ecodive' = cleanup(ecodive::variance_adjusted_unifrac(hmp50_mtx, hmp50_tree)) )
})
)
unifrac_res <- bind_rows(
mutate(u_unifrac_res, `UniFrac Variant` = 'Unweighted'),
mutate(w_unifrac_res, `UniFrac Variant` = 'Weighted'),
mutate(wn_unifrac_res, `UniFrac Variant` = 'Weighted Normalized'),
mutate(g_unifrac_res, `UniFrac Variant` = 'Generalized'),
mutate(va_unifrac_res, `UniFrac Variant` = 'Variance Adjusted') ) %>%
mutate(Package = as.character(expression)) %>%
select(Package, `UniFrac Variant`, median, mem_alloc) %>%
arrange(Package)
print(unifrac_res, n = 21)
#> # A tibble: 21 × 4
#> Package `UniFrac Variant` median mem_alloc
#> <chr> <chr> <bch:tm> <bch:byt>
#> 1 GUniFrac Unweighted 77.43ms 113.4MB
#> 2 GUniFrac Weighted Normalized 78.29ms 92.15MB
#> 3 GUniFrac Generalized 75.25ms 92.15MB
#> 4 abdiv Unweighted 14.36s 20.05GB
#> 5 abdiv Weighted 13.61s 20.02GB
#> 6 abdiv Weighted Normalized 13.59s 20.03GB
#> 7 abdiv Generalized 13.72s 20.18GB
#> 8 abdiv Variance Adjusted 15.57s 24.46GB
#> 9 ampvis2 Unweighted 3.31s 53.44MB
#> 10 ampvis2 Weighted 3.25s 49.2MB
#> 11 ampvis2 Weighted Normalized 3.3s 49.34MB
#> 12 ecodive Unweighted 3.94ms 1.01MB
#> 13 ecodive Weighted 3.46ms 779.72KB
#> 14 ecodive Weighted Normalized 3.82ms 779.73KB
#> 15 ecodive Generalized 5.59ms 799.3KB
#> 16 ecodive Variance Adjusted 4.12ms 779.73KB
#> 17 phyloregion Unweighted 6.35ms 146.44MB
#> 18 phyloseq Unweighted 290.5ms 50.87MB
#> 19 phyloseq Weighted 287.27ms 49.4MB
#> 20 phyloseq Weighted Normalized 291.59ms 49.55MB
#> 21 picante Unweighted 2.43s 1.79GB
# How much faster and more memory efficient is ecodive?
plyr::ddply(unifrac_res, as.name('UniFrac Variant'), function (x) {
x[['median']] <- as.numeric(x[['median']])
x[['mem_alloc']] <- as.numeric(x[['mem_alloc']])
ecodive <- as.list(x[x[['Package']] == 'ecodive',])
x <- x[x[['Package']] != 'ecodive',]
data.frame(
speed = paste0(paste(collapse=' - ', round(range(x$median / ecodive$median))), 'x'),
memory = paste0(paste(collapse=' - ', round(range(x$mem_alloc / ecodive$mem_alloc))), 'x') )
})
#> UniFrac Variant speed memory
#> 1 Generalized 13 - 2454x 118 - 26475x
#> 2 Unweighted 2 - 3647x 50 - 20248x
#> 3 Variance Adjusted 3776 - 3776x 32891 - 32891x
#> 4 Weighted 83 - 3931x 65 - 26922x
#> 5 Weighted Normalized 20 - 3555x 65 - 26934x
ggplot(unifrac_res, aes(x = median, y = mem_alloc)) +
geom_point(aes(shape = `UniFrac Variant`), size = 2) +
geom_label_repel(
data = ~subset(., `UniFrac Variant` == 'Unweighted'),
mapping = aes(label = Package),
box.padding = 0.4,
min.segment.length = Inf ) +
scale_shape(solid = FALSE) +
scale_y_log10(labels = scales::label_bytes()) +
labs(
title = 'UniFrac Implementations',
subtitle = 'All-vs-all 50 sample benchmarking on six CPU cores',
x = 'Median Calculation Time (log scale; n=10)',
y = 'Memory Allocated\n(log scale; n=1)' ) +
theme_bw(base_size = 12) +
theme(axis.title.x = element_text(margin = margin(t = 10)))
ecodive demonstrates substantial performance gains for
UniFrac, being 2 to 3,900x faster and using 50 - 32,000x less memory.
