--- title: "Output Consistency and Computational Cost Relative to Existing Tools" author: "Junhui Li, Lihua Julie Zhu" date: "`r Sys.Date()`" output: rmarkdown::html_vignette vignette: > %\VignetteIndexEntry{Output Consistency and Computational Cost Relative to Existing Tools} %\VignetteEngine{knitr::rmarkdown} %\VignetteEncoding{UTF-8} --- ```{r setup, include=FALSE} knitr::opts_chunk$set( echo = TRUE, message = FALSE, warning = FALSE, fig.align = "center", # html_vignette inherits fig.retina = 2 from rmarkdown, which renders every # figure at twice the nominal resolution and then scales it down in the # browser. The pixel count, and therefore the size of the base64 blob # embedded in the self-contained HTML, is four times larger for no visible # gain in a vignette. Setting it to 1 is what keeps the installed size of # doc/ within the limit R CMD check reports on. fig.retina = 1, dpi = 72, fig.width = 9, fig.height = 4.2 ) ``` # What is compared, and what is not Three implementations of the nearest-neighbor model are compared on identical input: **TmCalculator**, **rmelting** (an R interface to the MELTING 5 Java engine, whose JAR files are supplied by `melting5jars`, so the two packages share one computational engine), and **`Bio.SeqUtils.MeltingTemp`** from Biopython. Two things are deliberately kept apart. **Output consistency** asks whether the three tools, given the same parameter set and the same reaction conditions, return the same melting temperature. This part is independent of language, hardware and intended use, and it is the part worth trusting most. **Computational cost** asks how long each takes and how much memory it holds. Here the comparison is of *three tools as distributed*, not of three languages: TmCalculator's nearest-neighbor loop is compiled C++, Biopython's `Tm_NN` is pure Python, and MELTING 5 is Java. That difference is an implementation choice by each project, and it is what a user experiences, but it should not be read as a statement about R, Python or Java. The comparison is restricted to a plain list of oligonucleotides, which is the one task all three were designed to perform. Benchmarking MELTING 5 or Biopython on millions of genomic windows would compare tools built for different purposes. The coordinate model that only TmCalculator provides is a capability difference and appears below as a feature table, not as a timing. # Aligning the three tools A speed comparison between tools computing different models is meaningless, so the model is pinned in all three. | | TmCalculator | Biopython | MELTING 5 | |---|---|---|---| | Nearest-neighbor set | `nn_table = "DNA_NN_SantaLucia_2004"` | `nn_table = mt.DNA_NN4` | `method.nn = "san04"` | | Salt correction | `salt_method = "SantaLucia1996"` | `saltcorr = 3` | `correction.ion = "san96"` | | Sodium | `Na = 50` (mM) | `Na = 50` (mM) | `Na.conc = 0.05` (M) | | Strand concentration | `dnac_high = 25`, `dnac_low = 25` (nM) | `dnac1 = 25`, `dnac2 = 25` (nM) | `nucleic.acid.conc` (M) | Two points about that table are worth stating explicitly, because both were found by testing rather than by reading documentation. The salt correction is SantaLucia et al. (1996), `12.5 * log10[Na+]`, rather than the `11.7 * log10[Na+]` form used elsewhere in this package. It is the only plain additive logarithmic correction available in all three tools; MELTING 5's `correction.ion` has no entry for the 11.7 form, so that form cannot serve as common ground however natural it is otherwise. MELTING 5 defaults to `method.nn = "all97"` (Allawi & SantaLucia 1997) and `correction.ion = "ahs01"` (von Ahsen 2001). Neither is what the other two tools are running. Left unset, they produced a systematic offset of about 2.4 °C that had nothing to do with the implementations. The `nucleic.acid.conc` convention is the one parameter the R interface does not document, so `inst/scripts/bench_crosstool.R --grid` sweeps the plausible conventions against the candidate ion corrections and reports which pair agrees. Choosing the combination with the smallest deviation is *not* the right criterion, and the script prints enough to see why: switching MELTING's ion correction from `san96` to `schlif` shifts the mean deviation by exactly `(16.6 - 12.5) * log10(0.05) = -5.334` °C at every concentration, which is how one confirms that `san96` really is the same formula rather than merely the best-scoring one. # Reproducing the measurements These chunks are not evaluated when the vignette is built: they need a Java runtime, `rmelting`, and a Python interpreter with Biopython, none of which is a dependency of this package. ```{r reproduce, eval=FALSE} ## 1. Confirm the three tools agree before timing anything. ## If they do not, the benchmark measures default-value differences. system("Rscript inst/scripts/bench_crosstool.R --calibrate --outdir bench200") ## 2. rmelting accepts one sequence per call, so it is given its own, ## smaller set of sizes; results merge into the same CSV. system(paste("Rscript inst/scripts/bench_crosstool.R", "--tools rmelting --sizes 100,1000 --reps 3 --outdir bench200")) system(paste("Rscript inst/scripts/bench_crosstool.R", "--tools TmCalculator,Biopython", "--sizes 1000,10000,100000 --reps 3 --outdir bench200")) ## 3. Summary table and figure. system("Rscript inst/scripts/plot_crosstool.R --indir bench200") ``` Each tool runs in its own process under GNU `time`. That is the only way to attribute peak resident set size to a tool rather than to the accumulated state of one long session, and it means a tool that hangs does not take the others with it. Sequences are 200 bp windows drawn from the *E. coli* K-12 MG1655 chromosome under a fixed seed, matching the window size used in the genome-wide case study. # Results ```{r load} ## 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)) > 0L ``` ```{r nodata, eval=!has_data, echo=FALSE, results='asis'} ## message() is swallowed by the global message = FALSE, which would leave a ## build with no data looking simply empty. Write the explanation into the ## document itself instead. cat("> **The benchmark results are not present in this build.**\n>\n", "> `system.file(\"extdata\", \"crosstool_bench.csv\")` returned nothing,\n", "> so the tables and figures below are omitted. To populate them, run\n", "> `inst/scripts/bench_crosstool.R`, copy `bench_crosstool.csv` and\n", "> `consistency.csv` into `inst/extdata/` as `crosstool_bench.csv` and\n", "> `crosstool_consistency.csv`, and **reinstall the package** --\n", "> `system.file()` resolves against the installed copy, not the source\n", "> tree.\n", sep = "") ``` ```{r load2, eval=has_data} S <- 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.") ``` ## Consistency ```{r consistency, eval=has_data} 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.") } ``` 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 ```{r figure, eval=has_data, fig.cap="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."} 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) } par(op) ``` ```{r decompose, eval=has_data} ## 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.") ## 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))) } ``` 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. ```{r walltime, eval=has_data, fig.width=7.5, fig.height=4.6, fig.cap="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."} 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) ## 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 = ","))) } par(op) ``` 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. # What the timings do not show `rmelting::melting()` accepts one sequence per call: ```{r batch, eval=FALSE} rmelting::melting(sequence = c("ACGTACGTACGTACGTACGTACGT", "GGCCGGCCGGCCGGCCGGCCGGCC"), nucleic.acid.conc = 1.25e-8, hybridisation.type = "dnadna", Na.conc = 0.05, method.nn = "san04", correction.ion = "san96") #> Error: 'sequence' should be a character vector of length 1. ``` Processing *n* sequences therefore requires *n* separate invocations of the Java engine. The gap between MELTING 5 and the other two tools reflects interface design at least as much as the cost of the arithmetic. It also does not scale linearly. Between one hundred and one thousand sequences the input grew tenfold and the time grew about fiftyfold, so the rate fell from roughly nineteen sequences per second to four. Extrapolating from a small run therefore *understates* the cost of a large one, and the tool is omitted from the time figures for two reasons rather than one: at the sizes the other tools were measured at, it is orders of magnitude slower, and two points that do not lie on a straight line cannot honestly be drawn through sizes at which nothing was measured. Its measured values remain in the table above. The capability differences behind that observation are not timings and belong in a table of their own. | | Genomic coordinate model | Batch input | Parallel interface | `GRanges` output | |---|---|---|---|---| | TmCalculator | yes | yes | region-level, any backend | yes | | rmelting (MELTING 5) | no | no | no | no | | `Bio.SeqUtils.MeltingTemp` | no | yes | no | no | # Limitations All measurements come from a single machine and are not a portability claim. Repeated runs of the same configuration have differed by nearly 30% on this hardware, which is why the range over repetitions is shown rather than a standard deviation, and why conclusions are drawn from ratios spanning orders of magnitude rather than from small differences. Peak resident set size is measured for the whole process. R, Java and Python reserve memory differently, so these values indicate the practical footprint of each workflow rather than the memory the calculation itself requires. The consistency comparison covers perfectly matched DNA duplexes under one parameter set and one salt correction. It does not cover mismatches, dangling ends, RNA or RNA/DNA hybrids, for which the three tools do not offer the same parameter sets. # Session Information ```{r session-info} sessionInfo() ```