The pipeline and its serial cost

chr1 is 248,956,422 bp. Non-overlapping 200 bp windows after N-trimming give 1,244,682 windows.

chr_len <- GenomeInfoDb::seqlengths(genome)[["chr1"]]

t_coord <- system.time({
  bins <- make_genomiccoord(bsgenome = pkg, chromosomes = "chr1",
                            window = 200L, slide = 200L,
                            start = 1, end = chr_len, strand = "+")
})

t_extract <- system.time({
  gr_batch <- to_genomic_ranges_fast(list(pkg_name = pkg, seq = bins),
                                     method = "preload_chr")
})

base::gc(reset = TRUE)
t_tm <- system.time({
  tm_chr1 <- tm_calculate(gr_batch, method = "tm_nn",
                          nn_table = "DNA_NN_SantaLucia_2004", Na = 50)
})
base::gc()   # "max used" = peak R memory during the Tm step

rbind(t_coord, t_extract, t_tm)[, "elapsed"]
## t_coord t_extract      t_tm
##     5.5      15.4      30.8     (freshly booted, otherwise idle machine)

Measured on the test machine (R --vanilla, freshly booted, no other applications):

On a machine with normal background load the same steps run 20-35% slower (we measured 7.3 / 19.3 / 41.5 s in an ordinary desktop session); comparisons within one section below were always measured in a single session.

Memory: base::gc() reported a peak of ~1.9 GB of R-managed heap during the Tm step, and a dedicated worker process running the whole chr1 pipeline reaches ~3.5-4 GB resident (sequences, complements, a character copy handed to the compiled core, GRanges metadata, plus loaded packages). A useful planning figure: budget ~4 GB per chr1-sized pipeline, tapering to well under 2 GB for small chromosomes.

Two extraction details matter at this scale. Pass the input as list(pkg_name = ..., seq = ...) so the genome package is loaded once, and use method = "preload_chr" for dense tiling of few chromosomes (getSeq()-based "vectorized" is better for windows scattered across many chromosomes).

Why parallelism inside a call does not accelerate a single chromosome

Versions before 1.1.0 accepted a BPPARAM argument in tm_calculate() that split the sequences into one chunk per worker and ran the Tm loop on a BiocParallel backend. It did exactly what it advertised, yet on chr1 it never helped, and the argument has been removed. The measurement that decided it (taken with the 1.0.x series; the call no longer runs):

system.time({
  tm_serial <- tm_calculate(gr_batch, method = "tm_nn", Na = 50)
})
## elapsed ~ 30 s

## 1.0.x only: BPPARAM no longer exists
system.time({
  tm_snow <- tm_calculate(gr_batch, method = "tm_nn", Na = 50,
                          BPPARAM = SnowParam(workers = 5))
})
## elapsed ~ 55 s on the idle machine (33-80 s across repeated sessions)
## -- NEVER faster than serial

The serial time is stable across sessions; the parallel time is not, because its added costs (spawning workers, serializing half a gigabyte of sequence, loading packages in each worker) are highly sensitive to machine state. The decomposition explains why even its best case cannot win. Of the ~30 s serial runtime, only the per-window Tm loop (~9 s with the compiled core) is divided among workers. The rest is inherently serial or is new cost created by parallelization:

This is Amdahl’s law with an unusually small parallel fraction (~30%): even with infinitely many workers the floor is ~29 s — essentially the serial runtime — and on a loaded machine the added costs balloon well past it. With the pure-R core of previous releases the loop dominated (minutes for an input of this size) and the same SnowParam(5) call gave a 2.8x speedup; compiling the loop moved the bottleneck into the serial part and made within-call parallelism unprofitable. With the compiled core, a plain serial call is the fastest way to run inputs up to roughly a million windows, which is why tm_calculate() now runs serially and leaves parallelism to the caller, at the level of regions (Part 2).

A stronger warning applies to MulticoreParam() on large inputs, and it carries over to the region-level parallelism of Part 2. Forked workers share the manager’s memory copy-on-write, but R’s garbage collector writes to every object header it marks, so each worker’s GC forces the kernel to physically duplicate inherited pages. On a chr1-scale input (measured with the pure-R core, where the effect is easiest to see), MulticoreParam(5) ran slower than serial while burning ~18 CPU-minutes of system time on page duplication. When you parallelize, prefer SnowParam().

Parallelize at the chromosome level

The serial preprocessing that defeats within-call parallelism is itself embarrassingly parallel across chromosomes: window building, extraction, filtering and assembly for chr2 are independent of chr1. So for genome-scale work, give each worker one whole chromosome and run the entire pipeline — coordinates, extraction, tm_calculate() with the serial default — inside the worker:

chrs <- paste0("chr", c(1:22, "X", "Y"))
## Better: sort largest-first for load balance -- see "Choosing the
## worker count" below for the one-liner.
n_workers <- 5   # see "Choosing the worker count" below

runtime <- system.time({
  res_list <- bplapply(chrs, function(chr, pkg) {
    ## SnowParam workers are fresh R processes: load packages HERE,
    ## not in the manager session.
    suppressPackageStartupMessages(library(TmCalculator))
    suppressPackageStartupMessages(library(pkg, character.only = TRUE))

    genome  <- get(pkg, envir = asNamespace(pkg))
    chr_len <- GenomeInfoDb::seqlengths(genome)[[chr]]

    bins <- make_genomiccoord(bsgenome = pkg, chromosomes = chr,
                              window = 200L, slide = 200L,
                              start = 1, end = chr_len, strand = "+",
                              verbose = FALSE)
    gr  <- to_genomic_ranges_fast(list(pkg_name = pkg, seq = bins),
                                  method = "preload_chr")
    out <- tm_calculate(gr, method = "tm_nn",
                        nn_table = "DNA_NN_SantaLucia_2004",
                        Na = 50)$gr        # serial inside the worker

    ## Drop sequence columns before returning: Tm/GC are what we keep,
    ## and this cuts per-chromosome serialization from ~500 MB to a few MB.
    out$sequence   <- NULL
    out$complement <- NULL

    ## Record this worker's peak resident memory (GB) for budgeting.
    attr(out, "worker_rss_gb") <- ps::ps_memory_info()[["rss"]] / 1e9
    out
  }, pkg = pkg, BPPARAM = SnowParam(workers = n_workers))

  tm_genome <- unlist(GenomicRanges::GRangesList(res_list))
})

runtime
sapply(res_list, attr, "worker_rss_gb")   # per-worker memory check
length(tm_genome)                         # 14,687,330 windows
summary(tm_genome$Tm)
## Min. 47.2  1st Qu. 68.9  Median 71.5  Mean 72.0  3rd Qu. 74.8  Max. 101.1

Unlike the old within-call BPPARAM, this scheme parallelizes everything: window building, extraction, filtering and assembly, not just the Tm loop. Per-task communication is a chromosome name in and a slim GRanges out, so Amdahl’s law has little serial residue to bite on.

The chromosome is not the only unit of work available, though, and the choice of unit turns out to matter more than the worker count. The next section defines three schemes and the one after it measures all three over the same worker range.

Three ways to divide the work

The three differ in how many tasks there are and in when a task is assigned to a worker.

Static dispatch is SnowParam(workers = n) left alone. BiocParallel pre-splits the task list into one contiguous chunk per worker before anything runs, so a worker that draws a light chunk finishes early and then sits idle for the rest of the run.

Dynamic dispatch is the same task list with tasks = length(X). The manager hands out one element at a time, so a worker that finishes immediately takes the next. Sorting the chromosomes largest first makes this effective: the long tasks are placed while there is still short work left to fill the ragged end of the run.

Segment dispatch keeps dynamic assignment and changes the unit. Instead of 24 whole chromosomes the genome is cut into 50 Mb intervals, giving 73 tasks, and the longest single task falls from 55 s to 12 s. No schedule can finish sooner than its longest task, so this lowers the floor that the other two schemes are pressed against.

sl   <- GenomeInfoDb::seqlengths(genome)[paste0("chr", c(1:22, "X", "Y"))]
chrs <- names(sort(sl, decreasing = TRUE))          # largest first

tasks_chrom <- lapply(chrs, function(ch)
  list(chr = ch, start = 1L, end = as.integer(sl[[ch]])))

## Segment boundaries must be multiples of `slide`, or the window grid
## shifts between segments and the result stops matching a whole-
## chromosome run. 50 Mb = 250,000 x 200 bp.
seg_size  <- 50e6
tasks_seg <- unlist(lapply(names(sl), function(ch) {
  st <- seq(1, sl[[ch]], by = seg_size)
  lapply(st, function(s)
    list(chr = ch, start = as.integer(s),
         end = as.integer(min(s + seg_size - 1, sl[[ch]]))))
}), recursive = FALSE)

length(tasks_chrom)   # 24
length(tasks_seg)     # 73

n <- 5
BPPARAM_static  <- SnowParam(workers = n)
BPPARAM_dynamic <- SnowParam(workers = n, tasks = length(tasks_chrom))
BPPARAM_segment <- SnowParam(workers = n, tasks = length(tasks_seg))

One correctness detail: a segment task must pass trim_N = "none" to make_genomiccoord(), while a whole-chromosome task should trim. Trimming inside a segment would move the window grid relative to the chromosome and the two schemes would no longer compute the same windows. Interior N windows are dropped afterwards in both cases, so the totals agree exactly.

The measured sweep

All 24 chromosomes, 14,687,330 windows, 200 bp non-overlapping, freshly booted, R --vanilla, three repetitions per configuration in each environment. Wall time is elapsed time including worker start-up, which is also reported on its own; workers start concurrently, so start-up contributes its longest instance rather than the sum. Speedup is the one-worker wall time of the same strategy and environment divided by the wall time of the configuration, and efficiency is that speedup divided by the worker count.

# The two sweeps ship with the package as the CSV files the benchmark
# scripts wrote, and everything below is derived from them. Nothing is
# transcribed, so the tables, the figure and the numbers quoted in the text
# cannot drift apart, and re-running a sweep means replacing a file.
cols <- c("env", "strategy", "n_workers", "n_tasks", "rep", "wall_s",
          "startup_s", "wall_compute_s", "work_s", "hard_floor", "max_rss_gb")
read_sweep <- function(file, env) {
  d <- utils::read.csv(system.file("extdata", file, package = "TmCalculator"),
                       stringsAsFactors = FALSE)
  stopifnot(length(unique(d$n_windows)) == 1L)   # identical windows throughout
  d$env <- env
  d[, cols]
}
raw <- rbind(read_sweep("bench_parallel_strategy.csv", "laptop"),
             read_sweep("bench_parallel_cluster.csv",  "node"))

# Median over repetitions; the observed range is kept for the wall clock,
# which is the quantity the range bars in the figure show.
sweep <- do.call(rbind, lapply(
  split(raw, list(raw$env, raw$strategy, raw$n_workers), drop = TRUE),
  function(d) data.frame(
    env = d$env[1], strategy = d$strategy[1], workers = d$n_workers[1],
    tasks = d$n_tasks[1], n_rep = nrow(d),
    wall_s = median(d$wall_s), wall_lo = min(d$wall_s), wall_hi = max(d$wall_s),
    startup_s = median(d$startup_s), compute_s = median(d$wall_compute_s),
    work_s = median(d$work_s), longest_task_s = median(d$hard_floor),
    rss_gb = median(d$max_rss_gb), stringsAsFactors = FALSE)))

# Derived rather than transcribed. Speedup is the one-worker wall clock of
# the SAME strategy in the SAME environment over the configuration's wall
# clock: segmenting changes how much total work there is, so a ratio formed
# against a configuration's own summed task times would credit a strategy
# for doing more work and the three could not be compared.
key    <- paste(sweep$env, sweep$strategy)
ser    <- sweep[sweep$workers == 1L, ]
serial      <- stats::setNames(ser$wall_s, paste(ser$env, ser$strategy))
serial_work <- stats::setNames(ser$work_s, paste(ser$env, ser$strategy))
sweep$speedup    <- serial[key] / sweep$wall_s
sweep$efficiency <- sweep$speedup / sweep$workers
sweep$work_ratio <- sweep$work_s / serial_work[key]
sweep$strategy <- factor(sweep$strategy, levels = c("static", "dynamic", "segment"))
sweep$env      <- factor(sweep$env, levels = c("laptop", "node"))
sweep <- sweep[order(sweep$env, sweep$strategy, sweep$workers), ]
rownames(sweep) <- NULL

sweep_table <- function(s, caption) {
  tab <- data.frame(
    Strategy   = ifelse(duplicated(s$strategy), "", as.character(s$strategy)),
    Workers    = s$workers,
    Tasks      = s$tasks,
    `Wall (s)` = sprintf("%.1f [%.1f-%.1f]", s$wall_s, s$wall_lo, s$wall_hi),
    `Start-up (s)` = round(s$startup_s, 1),
    `Total task time (s)` = round(s$work_s, 1),
    `Work ratio` = round(s$work_ratio, 2),
    `Longest task (s)` = round(s$longest_task_s, 1),
    Speedup    = round(s$speedup, 2),
    Efficiency = round(s$efficiency, 2),
    `Peak RSS per worker (GB)` = round(s$rss_gb, 2),
    check.names = FALSE, stringsAsFactors = FALSE)
  knitr::kable(tab, row.names = FALSE, caption = caption)
}
sweep_table(sweep[sweep$env == "laptop", ],
  paste("Laptop (6-core Intel Core i7, 16 GB): task-partitioning strategy",
        "and worker count. Wall time is the median of three repetitions with",
        "the observed range in brackets."))

sweep_table(sweep[sweep$env == "node", ],
  paste("Compute node (two 20-core Intel Xeon Gold 6230, six slots, 48 GB",
        "limit, shared with other jobs): the same sweep."))

# Three panels, one argument, now made twice. Panel A shows that on the
# laptop every strategy turns over before the cores run out, which invites
# the usual explanation of a ragged schedule. Panel B rules that out: the
# total task time grows with the worker count, so the later configurations
# are doing more work, not dividing a fixed amount badly. Panel C gives the
# reason and the remedy together. The node, overlaid, shows what happens to
# each of the three when memory stops being scarce.
pal  <- c(static = "#1B5E9C", dynamic = "#C0392B", segment = "#5C6B73")
ltys <- c(laptop = "solid", node = "22")   # "22": a short dash, legible in a key
pchs <- c(laptop = 16, node = 1)
lv   <- levels(sweep$strategy)
ev   <- levels(sweep$env)
wk   <- sort(unique(sweep$workers))

op <- par(mfrow = c(1, 3), mar = c(4.4, 4.5, 2.2, 0.8), las = 1,
          cex = 0.8, mgp = c(2.8, 0.7, 0))

series <- function(col, ylab, ref = NULL, ref_lab = NULL, ymax = NULL,
                   diagonal = FALSE) {
  plot(NA, xlim = range(wk),
       ylim = c(0, if (is.null(ymax)) max(sweep[[col]]) * 1.05 else ymax),
       bty = "n", xaxt = "n", xlab = "Workers", ylab = ylab)
  axis(1, at = wk)
  if (diagonal) {
    abline(a = 0, b = 1, lty = 2, col = "grey55")
    text(max(wk), max(wk), "linear", adj = c(1.1, -0.4), cex = 0.75,
         col = "grey40")
  }
  if (!is.null(ref)) {
    abline(h = ref, lty = 2, col = "grey55")
    text(max(wk), ref, ref_lab, adj = c(1.1, -0.5), cex = 0.75, col = "grey40")
  }
  for (e in ev) for (st in lv) {
    d <- sweep[sweep$env == e & sweep$strategy == st, ]
    lines(d$workers, d[[col]], col = pal[st], lwd = 2, lty = ltys[e])
    points(d$workers, d[[col]], col = pal[st], pch = pchs[e], lwd = 1.6)
  }
}

series("speedup", "Speedup", ymax = max(wk) * 1.02, diagonal = TRUE)
mtext("A", side = 3, adj = 0, font = 2, line = 0.7, cex = 1.05)
# One legend, two columns: environment keys on the left, strategy colours
# on the right. legend() fills column-major, so the shorter column is padded
# with a blank entry whose lty, pch and col are NA.
legend("topleft", bty = "n", cex = 0.8, ncol = 2, seg.len = 2.6,
       legend = c(ev, "", lv),
       lty = c(unname(ltys[ev]), NA, rep("solid", length(lv))),
       pch = c(unname(pchs[ev]), NA, rep(NA, length(lv))),
       col = c("grey25", "grey25", NA, unname(pal[lv])), lwd = 2)

series("work_s", "Total task time (s)",
       ref = median(sweep$work_s[sweep$workers == 1]), ref_lab = "serial")
mtext("B", side = 3, adj = 0, font = 2, line = 0.7, cex = 1.05)

series("rss_gb", "Peak resident set size per worker (GB)")
mtext("C", side = 3, adj = 0, font = 2, line = 0.7, cex = 1.05)

Parallel performance under three task-partitioning strategies in two environments. Colour is the strategy; solid lines with filled points are the laptop, dashed lines with open points the compute node. (A) Speedup against worker count, relative to the one-worker run of the same strategy and environment; the dashed grey line marks linear speedup. (B) Total task time, summed over all tasks and measured inside the workers; values above the serial reference indicate that a configuration performed more work than the serial run rather than dividing the same work among more processes. (C) Peak resident set size of the heaviest worker.

par(op)

Three readings.

Dispatch order matters at every worker count. Static and dynamic run the identical 24 tasks; the only difference is that static commits them to workers up front. At four workers that costs 60 s, a fifth of the run.

Segmenting is worth more than dispatch order. On the laptop, segment dispatch is faster than dynamic at every worker count from three upwards, and its best time, 189 s, is 29 s better than dynamic’s best and 66 s better than static’s. The mechanism is in the “longest task” column: 55 s of chr1 cannot be divided by any scheduler, and at five workers that single task is nearly a third of the whole run.

Segmenting halves the memory. A worker holding a 50 Mb segment needs about 1.9 GB against 3.7 GB for one holding chr1, and peak memory is what decides how many workers a machine can run at all. On this 16 GB laptop segment dispatch is both the fastest scheme and the one that leaves the most headroom, which is an unusually convenient combination.

On the laptop, all three turn over well before the cores run out. Six physical cores are available, yet static peaks at 6, dynamic at 4 and segment at 5, and none reaches half of linear speedup. The next section explains why, the reason is not the one that first suggests itself, and the section after it shows the same sweep on a machine where the reason does not apply.

Why the speedup saturates

The obvious explanation for a curve that flattens is a ragged schedule: one long task holding up the end while other workers idle. The “total task time” column rules that out. It is the sum of the per-task times, measured inside the workers and excluding their package loading, so it is the amount of work the run actually performed. If the schedule were merely uneven, that column would stay flat and only the wall clock would suffer.

It does not stay flat. From one worker to six it grows by 53% under static dispatch, 88% under segment dispatch and 152% under dynamic dispatch. The later configurations are not dividing a fixed amount of work badly, they are doing more of it: the same chr1 task that takes 53 s alone takes 127 s with five siblings competing for memory bandwidth and last-level cache. Each worker is individually slower, and past a certain point adding one slows the others by more than it contributes.

This also explains the ordering. Dynamic dispatch inflates worst because its tasks are the largest resident objects, so the workers contend most for cache; segmenting reduces the resident footprint per worker and inflates less. The effect is a property of the machine rather than of the software, so a machine with more memory should show the turnover later, or not at all. That is a prediction, and the next section tests it. Rerun the sweep on your own hardware rather than transplanting these worker counts.

The work ratio column is the one to watch when you do. As long as it stays near 1.0, adding workers is buying real parallelism; once it passes about 1.3 the run is mostly paying for its own contention.

The same sweep on a compute node

The sweep was repeated with nothing changed but the machine: the same window grid, parameter set, worker counts and three repetitions, on a cluster node with two 20-core Intel Xeon Gold 6230 CPUs at 2.10 GHz, hyper-threading disabled, and 376 GB of memory, of which the job requested six slots and 48 GB. The node was shared with other jobs, so the run-to-run ranges are reported rather than suppressed; they are wide at one to three workers and narrow from four upwards.

Two things about that setup need saying. The job could not be given the laptop’s 16 GB. The scheduler enforces the sum of resident set size over every process in the job, and six socket workers plus their manager occupy about 23 GB, so a 16 GB request was killed after half an hour; 48 GB, three times the laptop and twice the measured peak, is what ran. And the workers share almost nothing with one another: summed proportional set size, which divides each shared page among its sharers, tracks summed RSS to within 0.2%, because each PSOCK worker attaches the packages and retrieves sequence for itself. The 23 GB is real, not an accounting artefact of shared mappings.

The result is the one the contention argument predicts, and it shows in every quantity that argument used:

The turnover is gone. Every curve on the node improves monotonically to six workers; dynamic dispatch reaches 128.9 s where on the laptop it was slowest at six of any configuration beyond three; and neither dynamic nor segment dispatch had saturated. The growth in total task time, which was the evidence that workers were slowing one another, falls to under a fifth. The longest task grows by a quarter instead of doubling. And the run-to-run spread, which on the laptop widened to 14 to 35 s at four to six workers for the two dynamic schemes, is 6 s or less on the node for the same configurations; timing that will not repeat is what paging looks like.

Segment dispatch is still the fastest, at 127.0 s, but its margin over dynamic dispatch falls from 29 s to 2 s. Bounding the working set is worth a great deal on a machine that is short of memory and almost nothing on one that is not. That is the practical form of the point made above: the ordering of the strategies is a property of the machine as much as of the partitioning.

One number runs the other way and is worth understanding rather than explaining away. Peak memory per worker is higher on the node, 3.83 GB against 2.15 GB for dynamic dispatch at six workers. R’s garbage collector reclaims on demand, so a process under no pressure keeps more allocated than one that is being squeezed; measured memory use is partly a function of memory available, and the per-worker figures in the tables are what the process settled at, not what it needed.

Where the extra work lands

The summary table reports only the longest task, which hides which task that was and whether the inflation is spread evenly. The benchmark also writes a per-task file, *_tasks.csv, with one row per task and repetition: chromosome, interval, window count, time, and the worker’s resident memory. The laptop’s is shipped with the package, compressed, and recovering the per-chromosome picture from it is a few lines:

tk <- utils::read.csv(system.file("extdata", "bench_parallel_strategy_tasks.csv.gz",
                                  package = "TmCalculator"))
tk <- tk[tk$strategy == "dynamic", ]        # one task per chromosome

## Median over repetitions, one column per worker count.
w <- stats::reshape(
  stats::aggregate(secs ~ chr + n_workers, data = tk, FUN = stats::median),
  idvar = "chr", timevar = "n_workers", direction = "wide")
names(w) <- sub("^secs\\.", "w", names(w))

w$inflation <- w$w6 / w$w1                  # 6 workers vs serial
w <- w[order(-w$w1), ]
knitr::kable(w, row.names = FALSE, digits = 1,
             caption = "Laptop, dynamic dispatch: median seconds per chromosome at each worker count.")

Read it two ways. Down the w1 column is the serial cost of each chromosome, which is what sets the floor no scheduler can go below: the largest single value is the hard limit on any chromosome-level schedule, and it is the number segmenting exists to reduce. Across a row is how much slower that same task became as workers were added, with inflation summarising it. A ratio near 1.0 means the task was indifferent to its neighbours; a large ratio means it was competing with them for memory bandwidth and cache.

The distinction matters for choosing a strategy. If inflation were uniform across chromosomes it would be a fixed tax and segmenting would not help with it. If instead the large chromosomes inflate most, then the same property that makes them stragglers also makes them the main source of contention, and cutting them into segments addresses both at once.

Segmenting one chromosome

Segment dispatch also answers the question Part 1 left open: how do you parallelize a single chromosome? Ship coordinates, not sequences. Each worker receives a (start, end) pair, a few dozen bytes, and runs the entire pipeline on its own segment. Every row of the Part 1 cost table becomes parallel except worker startup, because the preprocessing happens inside the workers and the half-gigabyte serialization never occurs.

chr_len <- GenomeInfoDb::seqlengths(genome)[["chr1"]]

## Segment length must be a multiple of `slide` so the window grid stays
## aligned across segment boundaries (50 Mb = 250,000 x 200).
seg_starts <- seq(1, chr_len, by = 50e6)
seg <- data.frame(start = seg_starts,
                  end   = pmin(seg_starts + 50e6 - 1, chr_len))

t_seg <- system.time({
  res_seg <- bplapply(seq_len(nrow(seg)), function(i, seg, pkg) {
    suppressPackageStartupMessages(library(TmCalculator))
    suppressPackageStartupMessages(library(pkg, character.only = TRUE))
    bins <- make_genomiccoord(bsgenome = pkg, chromosomes = "chr1",
                              window = 200L, slide = 200L,
                              start = seg$start[i], end = seg$end[i],
                              strand = "+", trim_N = "none",
                              verbose = FALSE)
    gr <- to_genomic_ranges_fast(list(pkg_name = pkg, seq = bins),
                                 method = "preload_chr")
    out <- tm_calculate(gr, method = "tm_nn",
                        nn_table = "DNA_NN_SantaLucia_2004", Na = 50)$gr
    out$sequence <- NULL; out$complement <- NULL
    out
  }, seg = seg, pkg = pkg, BPPARAM = SnowParam(workers = 5))

  tm_chr1_seg <- sort(unlist(GenomicRanges::GRangesList(res_seg)))
})

t_seg["elapsed"]      # 38 s measured, vs ~52 s for the serial FULL pipeline
length(tm_chr1_seg)   # 1,152,300

Mind the comparison scope: 38 s covers the entire pipeline (windows, extraction, Tm), so the serial reference is the ~52 s pipeline total from Part 1 — a 1.35x speedup — not the 30 s Tm-only figure. Modest, but it is the only single-chromosome parallelization that wins at all (within-call SnowParam lost to serial above), and the result is verifiably identical: 1,152,300 windows equals the serial run’s 1,244,682 generated minus its 92,382 N-gap windows, exactly. That equivalence is by construction — trim_N = "none" keeps the window grid identical across segments (the “Skipped … regions” warning from tm_calculate() is expected; hg38 chr1 carries an ~18 Mb interior centromeric gap), and the 50 Mb segment length is a multiple of slide.

The gain is capped well below worker-count because each worker pays fixed costs the serial run pays once: loading packages and reading the same 249 Mb chromosome to extract its segment, with five workers contending for the same file. Per-worker memory stays moderate, about 1 to 2 GB, since only a 50 Mb segment is windowed at a time.

Choosing the strategy and the worker count

Take segment dispatch. It was fastest at every worker count from three upwards, it halves the per-worker memory, and the memory saving is what lets a machine run more workers in the first place. Chromosome-level dispatch is worth keeping only if the downstream analysis wants whole chromosomes back as single objects.

For the worker count, two constraints, take the minimum:

mem_gb    <- 16                                     # your machine
cores     <- parallel::detectCores(logical = FALSE) # physical cores
per_worker_gb <- 2                                  # 50 Mb segment task
n_workers <- min(cores - 1L, floor((mem_gb - 4) / per_worker_gb))
n_workers
## 16 GB / 6 cores  -> 5 workers, which is what the sweep found fastest
## 32 GB / 8 cores  -> 7 workers
## 64 GB / 10 cores -> 9 workers

Use 2 GB per worker for segment tasks and 4 GB for whole-chromosome tasks; those are the measured peaks, not estimates. Do not count hyperthreads. The compiled core is CPU-bound and logical cores add little beyond scheduler overhead, and the laptop sweep shows the run is already contending for cache at five real cores. On the compute node the memory term is nowhere near binding, 22 workers at 48 GB, so the core term governs; with six slots the formula gives five, and the sweep was still improving at six. That is the case the formula is built for: it is a floor to start from, and the sweep is what sets the number.

Treat the memory bound as hard. When workers do collide the system swaps, and a swapping parallel run is slower than serial: we measured 19 GB of swap and a wedged machine in exactly this state, aggravated by leftover workers from interrupted runs.

Two smaller settings that the sweep shows are not optional. Sort the tasks largest-first, because the scheduler deals them out in order and a long task placed last has nothing left to overlap with. And pass tasks = length(X) to SnowParam(); leaving it at the default is what makes a run “static”, and that cost 60 s at four workers here for no benefit at all.

sl   <- GenomeInfoDb::seqlengths(genome)
chrs <- names(sort(sl[paste0("chr", c(1:22, "X", "Y"))], decreasing = TRUE))

Reproducing the sweep

The whole measurement is scripted, including the per-task instrumentation that separates worker start-up from compute and records each worker’s peak resident memory:

# laptop, or any single machine
Rscript inst/scripts/bench_parallel_strategy.R \
  --workers 1,2,3,4,5,6 --reps 3 --outfile bench_parallel_strategy.csv

# cluster node under LSF. Edit the conda environment, queue and CPU model
# in the #BSUB header first; the script activates the environment itself,
# checks that the package and hg38 are installed before it starts, and
# writes results/ under the directory it was submitted from.
bsub < inst/scripts/bench_parallel_cluster.lsf

# the two-environment table and figure of the manuscript
Rscript inst/scripts/make_figure5_two_env.R \
  --csv-hpc results/bench_parallel_cluster.csv \
  --csv-mac bench_parallel_strategy.csv

Include one worker in the sweep. It is the baseline every speedup is divided by, and BiocParallel runs a single worker inside the manager rather than spawning one, so that point measures the serial cost with nothing to serialise. Without it the only available ratios are against each strategy’s own summed task times, which are not comparable across strategies because segmenting changes how much total work there is.

For a cluster, inst/scripts/bench_parallel_cluster.R adds the three things a shared machine needs: it reads the allocated core count from the scheduler rather than from detectCores(), which on a shared node reports the whole machine and oversubscribes badly; it stages the genome package on node-local disk, because BSgenome reads the .2bit at random offsets and several workers hitting a shared filesystem measure storage latency instead of software; and it sweeps physical cores only. The submission file pins one CPU model, because a cluster whose nodes span several processor families would otherwise return a range that measures the scheduler; requests six slots and 48 GB, for the reasons given above; records the node’s load before and after the sweep; and samples summed RSS and PSS every ten seconds so the shared-memory question can be answered from the log rather than assumed.

Housekeeping

Four operational notes. First, benchmark in a session that did not restore a saved workspace ([Workspace loaded from ~/.RData] at startup means gigabytes of old objects are already resident); disable “Restore .RData at startup” in RStudio or start R --vanilla. Second, SnowParam() workers load the installed package: after devtools::load_all() development changes, run devtools::install() before benchmarking, and note that load_all() compiles C++ at -O0 — time only the installed build. The same applies to devtools::build_vignettes(), which evaluates vignette code against a -O0 debug build (timings come out ~6x slow); set options(pkg.build_extra_flags = FALSE) to disable the debug flags, or rely on R CMD build / devtools::check(), which install with normal optimization before running vignettes. Changing the flag option does not trigger recompilation — stale -O0 object files in src/ are reused as long as sources are unchanged — so run pkgbuild::clean_dll() after changing it. Third, if a parallel run errors or is interrupted, workers can linger and hold gigabytes; check for orphaned R processes (Activity Monitor on macOS, ps aux | grep R elsewhere) before rerunning. Fourth, worker message() output is buffered until each task finishes — a silent console does not mean a hung run; watch the worker processes’ CPU instead.