Results
## The measurements are shipped with the package rather than recomputed here,
## since reproducing them needs a Java runtime, rmelting and a Python
## interpreter with Biopython, none of which is a dependency. If the file is
## absent -- a source tree in which the benchmark has not been run -- the
## remaining sections are skipped rather than failing the build.
f <- system.file("extdata", "crosstool_bench.csv", package = "TmCalculator")
has_data <- nzchar(f) && file.exists(f) &&
length(readLines(f, n = 1L, warn = FALSE)) > 0LS <- utils::read.csv(f, stringsAsFactors = FALSE)
S <- S[S$ok & !is.na(S$compute_s), , drop = FALSE]
## The benchmark file accumulates across invocations, so it can hold rows
## taken before and after a change to the package. Averaging those together
## gives a table that describes no build at all, so refuse instead.
if ("pkg_version" %in% names(S)) {
vs <- sort(unique(stats::na.omit(S$pkg_version)))
if (length(vs) > 1L)
stop("crosstool_bench.csv mixes TmCalculator versions (",
paste(vs, collapse = ", "), "); re-run the benchmark with --fresh.")
} else {
stop("crosstool_bench.csv predates the split of start-up from file I/O. ",
"Re-run inst/scripts/bench_crosstool.R and copy the result into ",
"inst/extdata/, then reinstall.")
}
## Median and range over repetitions. With three repetitions a standard
## deviation is not a meaningful estimate, and repeated runs of the same
## configuration on a laptop have differed by nearly 30%, so the spread is
## reported as the range actually observed.
agg <- do.call(rbind, lapply(split(S, list(S$tool, S$n), drop = TRUE), function(d) {
data.frame(tool = d$tool[1], n = d$n[1], reps = nrow(d),
compute_med = stats::median(d$compute_s),
compute_lo = min(d$compute_s),
compute_hi = max(d$compute_s),
seqs_per_s = d$n[1] / stats::median(d$compute_s),
rss_med = stats::median(d$rss_gb),
rss_lo = min(d$rss_gb), rss_hi = max(d$rss_gb),
startup_med = stats::median(d$startup_s),
## Reading the input and writing the results exist only because
## each tool runs as a separate process, so that peak memory can
## be attributed to it. Reported, never added to start-up, never
## drawn: a user passes data in memory.
io_med = stats::median(d$io_s),
stringsAsFactors = FALSE)
}))
agg <- agg[order(agg$tool, agg$n), ]
## Start-up cannot depend on the input. If it does, something that scales
## with n is being counted as start-up.
for (t in unique(agg$tool)) {
d <- agg[agg$tool == t, ]
if (nrow(d) > 1L && max(d$startup_med) / min(d$startup_med) > 1.5)
warning(t, " start-up varies ",
sprintf("%.1fx", max(d$startup_med) / min(d$startup_med)),
" across input sizes; it should be constant.", call. = FALSE)
}
knitr::kable(agg, digits = 4, row.names = FALSE,
caption = "Median and range over repetitions at each input size.")| tool | n | reps | compute_med | compute_lo | compute_hi | seqs_per_s | rss_med | rss_lo | rss_hi | startup_med | io_med |
|---|---|---|---|---|---|---|---|---|---|---|---|
| Biopython | 1e+02 | 3 | 0.0210 | 0.0205 | 0.0213 | 4761.9048 | 0.0103 | 0.0103 | 0.0104 | 0.0941 | 0.0006 |
| Biopython | 3e+02 | 3 | 0.0625 | 0.0618 | 0.0628 | 4800.0000 | 0.0104 | 0.0104 | 0.0104 | 0.0958 | 0.0007 |
| Biopython | 1e+03 | 3 | 0.2136 | 0.2047 | 0.2137 | 4681.6479 | 0.0107 | 0.0106 | 0.0107 | 0.0960 | 0.0014 |
| Biopython | 3e+03 | 3 | 0.6207 | 0.6145 | 0.6228 | 4833.2528 | 0.0115 | 0.0115 | 0.0118 | 0.0990 | 0.0031 |
| Biopython | 1e+04 | 3 | 2.0701 | 2.0672 | 2.0781 | 4830.6845 | 0.0145 | 0.0145 | 0.0146 | 0.0995 | 0.0094 |
| Biopython | 3e+04 | 3 | 6.6377 | 6.2713 | 7.1749 | 4519.6378 | 0.0243 | 0.0241 | 0.0244 | 0.1529 | 0.0285 |
| Biopython | 1e+05 | 3 | 21.2312 | 20.6597 | 21.2901 | 4710.0494 | 0.0553 | 0.0549 | 0.0558 | 0.1322 | 0.0920 |
| Biopython | 3e+05 | 3 | 62.8029 | 61.6619 | 63.2765 | 4776.8495 | 0.1500 | 0.1499 | 0.1500 | 0.1683 | 0.2708 |
| TmCalculator | 1e+02 | 3 | 0.0960 | 0.0960 | 0.1230 | 1041.6667 | 0.3392 | 0.3392 | 0.3401 | 2.5820 | 0.0030 |
| TmCalculator | 3e+02 | 3 | 0.1050 | 0.1040 | 0.1330 | 2857.1429 | 0.3404 | 0.3401 | 0.3405 | 2.5240 | 0.0030 |
| TmCalculator | 1e+03 | 3 | 0.1330 | 0.1320 | 0.1400 | 7518.7970 | 0.3430 | 0.3426 | 0.3439 | 2.4360 | 0.0050 |
| TmCalculator | 3e+03 | 3 | 0.2140 | 0.2130 | 0.2150 | 14018.6916 | 0.3436 | 0.3434 | 0.3470 | 2.4240 | 0.0100 |
| TmCalculator | 1e+04 | 3 | 0.4940 | 0.4930 | 0.5040 | 20242.9150 | 0.3509 | 0.3509 | 0.3525 | 2.4310 | 0.0280 |
| TmCalculator | 3e+04 | 3 | 1.3260 | 1.3240 | 1.3290 | 22624.4344 | 0.3522 | 0.3518 | 0.3558 | 2.4170 | 0.1020 |
| TmCalculator | 1e+05 | 3 | 4.5290 | 4.5160 | 4.5470 | 22079.9293 | 0.3891 | 0.3850 | 0.3896 | 2.4680 | 0.2750 |
| TmCalculator | 3e+05 | 3 | 13.6470 | 13.3390 | 13.8150 | 21982.8534 | 0.5793 | 0.5785 | 0.5804 | 2.3420 | 1.3740 |
| rmelting | 1e+02 | 3 | 5.3290 | 5.3010 | 5.3730 | 18.7652 | 0.1810 | 0.1802 | 0.1822 | 0.5500 | 0.0020 |
| rmelting | 3e+02 | 3 | 28.4790 | 28.3620 | 28.5170 | 10.5341 | 0.2027 | 0.2020 | 0.2033 | 0.5380 | 0.0030 |
| rmelting | 1e+03 | 3 | 263.2020 | 254.5180 | 267.8610 | 3.7994 | 0.2824 | 0.2796 | 0.3149 | 0.7270 | 0.0060 |
Consistency
fc <- system.file("extdata", "crosstool_consistency.csv", package = "TmCalculator")
if (nzchar(fc)) {
knitr::kable(utils::read.csv(fc), digits = 10, row.names = FALSE,
caption = "Deviation in Tm relative to TmCalculator on identical input.")
}| tool | n | mean_dTm | max_abs_dTm | pearson_r |
|---|---|---|---|---|
| rmelting | 100 | 0.4787223867 | 1.350259800 | 0.9975431 |
| Biopython | 100 | 0.0000000001 | 0.000000005 | 1.0000000 |
Under an identical nearest-neighbor set and salt correction,
TmCalculator and Biopython return the same value to within
floating-point representation. MELTING 5 differs by a fraction of a
degree in a sequence-dependent way. Recovering enthalpy and entropy from
Tm measured at two strand concentrations
(bench_crosstool.R --thermo) localises that difference:
MELTING’s enthalpies are offset by a near-constant 0.30 kcal/mol, which
is a difference in the duplex initiation term rather than in the
stacking sum, while the entropies differ by an amount that is not
constant. Nominally identical parameter sets are therefore not
implemented identically across tools.
This comparison is also how an error in TmCalculator itself was
found: four of the six reverse-complement rows added to every
nearest-neighbor table had been transposed. Before that fix the
deviation from Biopython was sequence-dependent and up to 0.09 °C; after
it, the two agree exactly. See NEWS.md.
Cost
tools <- sort(unique(agg$tool))
pal <- c("#1B5E9C", "#C0392B", "#5C6B73")[seq_along(tools)]; names(pal) <- tools
pch <- c(16, 17, 15)[seq_along(tools)]; names(pch) <- tools
draw_range <- function(x, lo, hi, col) {
v <- is.finite(lo) & is.finite(hi) & hi / pmax(lo, 1e-12) > 1.02
if (any(v)) arrows(x[v], lo[v], x[v], hi[v], code = 3, angle = 90,
length = 0.03, col = col)
}
op <- par(mfrow = c(1, 2), mar = c(4.3, 4.4, 2.2, 0.8), las = 1, cex = 0.85)
xr <- range(agg$n)
## A tool measured at only the smallest sizes is left off the time axis: its
## cost is orders of magnitude larger, which would flatten the others onto the
## axis, and drawing two points as a line invites interpolation through sizes
## at which it was never run.
tt <- sort(unique(agg$tool[agg$n == max(agg$n)]))
aggt <- agg[agg$tool %in% tt, ]
plot(NA, xlim = c(0, max(aggt$n)), ylim = c(0, max(aggt$compute_hi) * 1.05),
bty = "n", xlab = "Sequences", ylab = "Compute time (s)", xaxt = "n")
axis(1, at = pretty(c(0, max(aggt$n))),
labels = format(pretty(c(0, max(aggt$n))), big.mark = ",",
scientific = FALSE, trim = TRUE))
mtext("A", side = 3, adj = 0, font = 2, line = 0.8, cex = 1.1)
for (t in tt) {
d <- aggt[aggt$tool == t, ]; d <- d[order(d$n), ]
lines(d$n, d$compute_med, col = pal[t], lwd = 2)
draw_range(d$n, d$compute_lo, d$compute_hi, pal[t])
points(d$n, d$compute_med, col = pal[t], pch = pch[t], cex = 1.05)
}
legend("topleft", bty = "n", legend = tt, col = pal[tt],
pch = pch[tt], lwd = 2, cex = 0.9)
## Expand to whole decades: range() alone puts the extreme points on the frame
yr2 <- range(c(agg$rss_lo, agg$rss_hi))
yr2 <- c(10^floor(log10(yr2[1])), 10^ceiling(log10(yr2[2])))
plot(NA, xlim = xr, ylim = yr2, log = "xy", bty = "n", yaxt = "n",
xlab = "Sequences", ylab = "Peak resident set size (GB)")
ticks <- 10^seq(log10(yr2[1]), log10(yr2[2]))
axis(2, at = ticks, labels = format(ticks, scientific = FALSE, drop0trailing = TRUE))
axis(2, at = as.numeric(outer(2:9, ticks)), labels = FALSE, tcl = -0.2)
mtext("B", side = 3, adj = 0, font = 2, line = 0.8, cex = 1.1)
for (t in tools) {
d <- agg[agg$tool == t, ]; d <- d[order(d$n), ]
if (nrow(d) > 1) lines(d$n, d$rss_med, col = pal[t], lwd = 2)
draw_range(d$n, d$rss_lo, d$rss_hi, pal[t])
points(d$n, d$rss_med, col = pal[t], pch = pch[t], cex = 1.1)
}
Compute time (A) and peak resident set size (B) against input size. Panel A is linear on both axes: above a few thousand sequences each tool is proportional to its input, so each is a straight line and the ratio of the slopes is the ratio of the throughputs, read directly off the picture. On logarithmic axes that ratio would carry the same visual weight as every other ratio in the plot, including the reversed one at the smallest input. Panel A omits any tool not measured across the whole range; panel B keeps all of them, the spread there being small enough to show.
## Two points at the top of the range separate the constant part of a call
## from the part that scales with the input.
fit <- do.call(rbind, lapply(split(agg, agg$tool), function(d) {
d <- d[order(d$n), ]
if (nrow(d) < 2L) return(NULL)
k <- nrow(d)
slope <- (d$compute_med[k] - d$compute_med[k - 1L]) / (d$n[k] - d$n[k - 1L])
data.frame(tool = d$tool[1],
fixed_cost_s = d$compute_med[k] - slope * d$n[k],
marginal_seqs_per_s = 1 / slope, stringsAsFactors = FALSE)
}))
knitr::kable(fit, digits = 3, row.names = FALSE,
caption = "Fixed cost per call and marginal throughput.")| tool | fixed_cost_s | marginal_seqs_per_s |
|---|---|---|
| Biopython | 0.445 | 4810.965 |
| TmCalculator | -0.030 | 21934.635 |
| rmelting | -72.117 | 2.982 |
## The size at which two tools take equal time follows from those two
## numbers, and is the one figure that keeps either end of panel A from
## being over-read.
if (nrow(fit) >= 2L) {
o <- order(fit$marginal_seqs_per_s, decreasing = TRUE)
a <- fit[o[1], ]; b <- fit[o[2], ]
crossover <- (a$fixed_cost_s - b$fixed_cost_s) /
(1 / b$marginal_seqs_per_s - 1 / a$marginal_seqs_per_s)
cat(sprintf("%s overtakes %s at about %s sequences.\n",
a$tool, b$tool, signif(crossover, 2)))
}## TmCalculator overtakes Biopython at about -2900 sequences.Panel A shows why a single input size is not enough. TmCalculator
pays a fixed cost per call – constructing the GRanges,
loading the parameter tables, assembling the result – that is unchanged
whether one sequence is submitted or two hundred thousand. Below the
crossover printed above, that cost dominates, and a benchmark run only
at that scale reports TmCalculator as the slower tool. Above it the two
lines separate in the other direction and settle at the marginal
throughputs in the table, which is what “how fast is the calculation”
actually means.
That ratio, however, is not what a user invoking a tool once from the command line experiences, because it excludes the cost of starting the process.
wall <- aggt[order(aggt$n, aggt$tool), ]
m <- rbind(compute = wall$compute_med, `start-up` = wall$startup_med)
colnames(m) <- paste(wall$tool, wall$n)
grp <- cumsum(c(1, diff(wall$n) != 0))
space <- ifelse(c(TRUE, diff(grp) != 0), 1.1, 0.18)
fill <- c("#1B5E9C", "#BFD3E6")
op <- par(mar = c(5.6, 4.6, 2.4, 0.8), las = 1, cex = 0.85)
bp <- barplot(m, space = space, col = fill, border = NA, ylab = "",
names.arg = rep("", ncol(m)), ylim = c(0, max(colSums(m)) * 1.18))
u <- par("usr")
text(bp, u[3] - 0.02 * (u[4] - u[3]), labels = wall$tool, srt = 40, adj = 1,
xpd = TRUE, cex = 0.72)
mtext(format(unique(wall$n), big.mark = ","), side = 1, line = 3.6,
at = tapply(bp, grp, mean), cex = 0.85)
mtext("Sequences", side = 1, line = 4.7, cex = 0.85)
mtext("Elapsed time (s)", side = 2, line = 3.2, las = 0, cex = 0.85)
legend("topleft", bty = "n", fill = fill, border = NA,
legend = c("compute", "start-up"), cex = 0.9)
Elapsed time for one invocation, separated into start-up and compute. Start-up is the constant part: an R session with the Bioconductor packages attached costs about two seconds whatever the input, which is more than an entire Biopython run at the smaller sizes. Only tools measured across the whole range are shown.
## Where the TOTAL, not the compute time, crosses over.
if (nrow(fit) >= 2L) {
o <- order(fit$marginal_seqs_per_s, decreasing = TRUE)
a <- fit[o[1], ]; b <- fit[o[2], ]
su <- vapply(list(a, b), function(z)
stats::median(agg$startup_med[agg$tool == z$tool]), numeric(1))
n_eq <- ((a$fixed_cost_s + su[1]) - (b$fixed_cost_s + su[2])) /
(1 / b$marginal_seqs_per_s - 1 / a$marginal_seqs_per_s)
cat(sprintf("Including start-up, %s overtakes %s at about %s sequences.\n",
a$tool, b$tool, format(signif(n_eq, 2), big.mark = ",")))
}## Including start-up, TmCalculator overtakes Biopython at about 11,000 sequences.The two crossovers differ by more than an order of magnitude, and a reader who takes one for the other will reach the wrong conclusion about which tool to use. Comparing only the calculation, TmCalculator overtakes Biopython at a few hundred sequences; comparing what one command actually takes, it does not overtake until roughly ten thousand, because about two seconds are spent attaching R and the Bioconductor packages before any sequence is read. Both numbers are correct, and which one applies depends on whether the process is started once and then works, as in a genome-wide run, or started afresh for every small job.
Panel B is the half of the comparison that does not favour a Bioconductor package. TmCalculator’s memory is essentially a constant floor: an R session with Biostrings and GenomicRanges attached starts near 0.8 GB whatever the input, and grows barely at all across three orders of magnitude. Biopython holds only the sequences and grows from almost nothing. At the sizes measured here Biopython is far lighter; the two would meet only at a much larger input. The same applies to start-up, which is about two seconds for R with the Bioconductor stack against a fraction of a second for a Python interpreter. Both are one-time costs, invisible in a genome-scale run and material in a loop over short sequences.