Introduction
tm_calculate() turns a genome, or any set of regions in
it, into a melting temperature profile: it tiles the regions into
windows, fetches each window’s sequence, computes a Tm, and returns the
lot as one object. The same call takes sequences you already have, so
there is one function to learn; on a mammalian genome it spreads the
work across processes.
This vignette answers three questions, in order:
- What can I give it? A BSgenome package, a FASTA
file, or sequences already in R, and any set of regions within
them.
- How does it divide the work? Into tasks, one region
or one segment each, dispatched to workers.
- How many workers should I use, and how long will it
take? Measured on the human genome in two hardware
environments.
Before you start. You need a BSgenome package for
the genome you want to profile; this vignette uses
BSgenome.Hsapiens.UCSC.hg38. BiocParallel
comes with TmCalculator, so there is nothing to install for the parallel
part; attach it when you want to name a backend.
BiocManager::install("BSgenome.Hsapiens.UCSC.hg38")
The hg38 chunks below are not evaluated when the vignette is built,
since they take minutes; the timing table and the figure are evaluated,
and are built from measurements that ship with the package.
Quick start
The smallest genome-scale call needs a genome and a window width:
library(TmCalculator)
hg38 <- "BSgenome.Hsapiens.UCSC.hg38"
tm <- tm_calculate(hg38, window = 200, slide = 200)$gr
tm
## GRanges object with 14687412 ranges and 2 metadata columns:
## seqnames ranges strand | Tm GC
tm_calculate() returns a TmCalculator
object: $gr is the profile and $options
records the model and tiling it was produced with.
That tiles every standard chromosome into non-overlapping 200 bp
windows and computes a nearest-neighbour Tm for each, in a single
process. The window width has no default at genome scale, and is not
meant to: it is the resolution of the profile, so it belongs in the call
rather than in a default that quietly decides how many rows come back.
Adding workers is one more argument:
tm <- tm_calculate(hg38, window = 200, slide = 200,
BPPARAM = BiocParallel::SnowParam(workers = 5))$gr
Same windows, same Tm values, about a third of the time. How the
genome is divided changes the runtime, not the result. See the sweep
below for how far adding workers gets you.
Everything else is a refinement of those two calls: which regions,
which thermodynamic model, and how many workers.
What can I give
it?
The sequence
source
The first argument says where sequence comes from. It is given by
name, not as a loaded object, because each worker opens the
source for itself.
tm_calculate(hg38, window = 200) # an installed BSgenome package
tm_calculate("contigs.fa.gz", window = 200) # a FASTA file, gzip is fine
tm_calculate(oligos) # a character vector of sequences
A vector of sequences is staged as a temporary FASTA under
tmpdir and deleted on exit, so that the workers read it
rather than receive it. That is worth doing for a large set, because it
moves window construction and result assembly into the workers as well;
for a handful of sequences the staging and the worker start-up cost more
than the calculation, and tm_calculate() is the right call.
On a cluster, set tmpdir to node-local scratch: the default
tempdir() is often a small partition.
The regions
regions says what to take from the source, in whichever
form is at hand.
## Chromosomes or records, by name or by number. On a BSgenome the chr
## prefix is added or removed to match the genome, so the same code works on
## UCSC and on Ensembl; FASTA record names are matched exactly.
tm_calculate(hg38, regions = 1:22, window = 200) # autosomes
tm_calculate(hg38, regions = paste0("chr", c(1:22, "X", "Y")), # no chrM
window = 200)
tm_calculate("contigs.fa.gz", regions = c("contig_7", "contig_9"),
window = 200)
## Coordinate intervals. Commas and scientific notation are accepted, so a
## number pasted out of a genome browser works as it stands.
tm_calculate(hg38, regions = c("chr1:1-10e6", "chrX:5,000,000-6,000,000"),
window = 200)
## A mixture, when some chromosomes are wanted whole and others in part.
tm_calculate(hg38, regions = c(1:20, "chr21:1-10e6", "X", "Y"), window = 200)
## A GRanges, when the regions come from an annotation.
prom <- promoters(genes(TxDb.Hsapiens.UCSC.hg38.knownGene),
upstream = 1000, downstream = 500)
tm_calculate(hg38, regions = prom, window = 50, slide = 25)
## A GRanges carrying its own sequences is a source rather than a query, and
## then regions selects by overlap: a seqname takes every range on it, an
## interval takes the ranges it meets. Whole ranges come back, not clipped
## pieces of them, since their sequences are already fixed.
tm_calculate(probes_gr, regions = "chr7")
tm_calculate(probes_gr, regions = "chr7:1-1e6")
## window = NULL, the default, gives one window per region, which is what
## short records call for: a FASTA of array probes, primers or synthetic
## oligos returns one Tm per record. It is refused for a region over 1 Mb,
## where a single Tm would mean nothing.
tm_calculate("probes.fa", window = NULL)
tm_calculate(oligos, BPPARAM = BiocParallel::SnowParam(5))
Three things are worth knowing before you trust the output.
The default includes chrM. With no regions, a
BSgenome source covers GenomeInfoDb::standardChromosomes(),
which for GRCh38 is the 24 assembled chromosomes and
the mitochondrion. Name the chromosomes explicitly if that does not
belong in your profile. The sweep later in this vignette uses
paste0("chr", c(1:22, "X", "Y")) for exactly that reason,
and its 14,687,330 windows are the 24 without chrM.
Whole chromosomes are trimmed, named regions are not. A
region that names a whole chromosome has its leading and trailing
assembly gaps trimmed, since a telomeric run of N carries no windows. A
region given by coordinate is tiled from the start you asked for.
Windows containing N are dropped either way.
Overlapping regions are not merged. They produce windows
that appear twice in the result, so tm_calculate() warns
rather than double-counting them quietly.
Array probes
Array manifests are not read directly, since their layout is vendor
and version specific. Take the probe sequences out with the package that
already understands the manifest, illuminaio,
minfi or sesame for Infinium arrays, then pass
them as a FASTA file or straight to tm_calculate().
How does it divide the
work?
A task is the unit one worker owns from start to
finish: it opens the source, builds its own windows, fetches its own
sequence and computes its own Tm. Only a name and a coordinate pair
cross between processes, so no sequence is ever serialized.
unit decides how tasks are cut.
tm_calculate(hg38, window = 200, unit = "segment", # the default: 73 tasks
segment_size = 50e6, BPPARAM = BiocParallel::SnowParam(5))
tm_calculate(hg38, window = 200, unit = "region", # 24 tasks
BPPARAM = BiocParallel::SnowParam(5))
The two calls return the same profile, row for row.
regions and window decide what is computed;
unit and segment_size decide only how that
work is handed out, and segment boundaries are held to multiples of
slide so the window grid cannot shift when the segmenting
changes.
Segments are the better default, for two reasons.
Chromosome 1 is an indivisible task of 249 Mb, so with one task per
chromosome the run cannot finish before chromosome 1 does, however many
workers are available. And a worker holding a 50 Mb segment needs less
memory than one holding a whole chromosome.
Prefer SnowParam() to
MulticoreParam(). 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 collection forces the
kernel to duplicate the inherited pages. On a genome-scale input that
has been measured running slower than a single process.
How many workers, and
how long?
The sweep below ships with the package: hg38 at 200 bp
non-overlapping windows, 14,687,330 windows, one to six workers, three
repetitions per configuration, on a six-core 16 GB laptop and on a
cluster compute node given six slots and the same 16 GB, so that the two
differ in their processors and not in their quota.
# Read the shipped summaries rather than transcribing numbers, so the table
# cannot drift from the measurements it describes.
read_sweep <- function(file, env) {
d <- utils::read.csv(system.file("extdata", file, package = "TmCalculator"),
stringsAsFactors = FALSE)
d$Environment <- env
d
}
sweep <- rbind(read_sweep("bench_hg38_laptop.csv", "Laptop"),
read_sweep("bench_hg38_cluster.csv", "Compute node"))
sweep <- sweep[order(sweep$Environment != "Laptop", sweep$n_workers), ]
knitr::kable(
data.frame(
Environment = sweep$Environment,
Workers = sweep$n_workers,
`Wall time (s)` = sprintf("%.1f [%.1f-%.1f]", sweep$wall_s, sweep$lo, sweep$hi),
Speedup = sprintf("%.2f", sweep$speedup),
`Peak RSS per worker (GB)` = sprintf("%.2f", sweep$peak_worker_gb),
check.names = FALSE),
row.names = FALSE,
caption = paste("Median of three repetitions, observed range in brackets.",
"Wall time includes worker start-up."))
Median of three repetitions, observed range in brackets. Wall
time includes worker start-up.
| Laptop |
1 |
590.3 [586.1-590.7] |
1.00 |
3.47 |
| Laptop |
2 |
336.8 [335.4-337.6] |
1.75 |
2.50 |
| Laptop |
3 |
239.6 [233.1-243.9] |
2.46 |
2.63 |
| Laptop |
4 |
197.4 [196.4-199.8] |
2.99 |
2.18 |
| Laptop |
5 |
179.7 [178.0-184.0] |
3.28 |
2.00 |
| Laptop |
6 |
179.8 [176.4-214.9] |
3.28 |
2.81 |
| Compute node |
1 |
494.5 [485.8-529.0] |
1.00 |
3.55 |
| Compute node |
2 |
277.4 [276.8-278.1] |
1.78 |
3.53 |
| Compute node |
3 |
197.0 [196.7-197.3] |
2.51 |
2.88 |
| Compute node |
4 |
154.3 [154.1-156.3] |
3.20 |
2.95 |
| Compute node |
5 |
129.9 [128.8-130.3] |
3.81 |
3.10 |
| Compute node |
6 |
114.4 [114.0-116.0] |
4.32 |
3.59 |
op <- par(mfrow = c(1, 2), mar = c(4.2, 4.4, 2.2, 0.8), mgp = c(2.6, 0.7, 0))
for (what in c("wall", "rss")) {
ys <- if (what == "wall") sweep$wall_s else sweep$peak_worker_gb
plot(range(sweep$n_workers), range(0, ys * 1.05), type = "n",
xlab = "Workers", xaxt = "n", adj = 0,
ylab = if (what == "wall") "Wall clock (s)" else "Peak resident size per worker (GB)",
main = if (what == "wall") "A" else "B")
axis(1, at = sort(unique(sweep$n_workers)))
for (e in unique(sweep$Environment)) {
d <- sweep[sweep$Environment == e, ]
y <- if (what == "wall") d$wall_s else d$peak_worker_gb
solid <- e == "Laptop"
if (what == "wall")
arrows(d$n_workers, d$lo, d$n_workers, d$hi, angle = 90, code = 3,
length = 0.03, col = "grey40")
lines(d$n_workers, y, lty = if (solid) 1 else 2, col = "grey20")
points(d$n_workers, y, pch = if (solid) 19 else 1, col = "grey20")
}
if (what == "wall")
legend("topright", c("Laptop, 6 cores", "Compute node, 6 slots"),
lty = c(1, 2), pch = c(19, 1), bty = "n", cex = 0.85, col = "grey20")
}
Wall-clock time and peak memory per worker
against worker count, on a six-core 16 GB laptop (solid, filled) and a
compute node given six slots and 16 GB (dashed, open). Points are
medians of three repetitions; bars give the observed range.
Two things in that table decide a worker count.
More workers stop helping, and where depends on
cores. The laptop is fastest at five workers, 179.7 s against
590.3 s in one process, and a sixth worker returns the same time with a
far wider spread across repetitions (176.4 to 214.9 s, against 10.9 s or
less at every smaller count). Five workers and the process that
dispatches them already occupy the laptop’s six cores, so the sixth has
none left. The node, whose six slots sit on a forty-core machine,
improves all the way to six workers and 114.4 s, its repetitions
differing by 2.2 s or less throughout.
Memory was not the limit on either machine. Peak
resident size per worker reached 3.5 GB in the one-worker run and stayed
between 2.0 and 3.6 GB thereafter, because no worker holds more than one
50 Mb task at a time. Both sweeps ran well inside 16 GB.
Fitting the wall times to Amdahl’s law gives
T(n) = 79.4 + 509.2/n on the laptop and
T(n) = 41.6 + 456.4/n on the node, both with R-squared of
at least 0.997. The divisible part is nearly the same on the two
machines, as it should be for the same work; what differs is the
constant, by a factor of 1.9, and the number of cores available to
divide among.
A rule that works on both machines:
workers = physical cores - 1
Confirm it on your own machine with a short sweep over one chromosome
rather than the whole genome. Note that worker start-up, roughly nine
seconds per call, is included in every timing, so a run that takes less
than a minute will understate the speedup:
for (n in 1:6) {
t <- system.time(
tm_calculate(hg38, regions = "chr21", window = 200, slide = 200,
method = "tm_nn", segment_size = 10e6,
BPPARAM = BiocParallel::SnowParam(n), verbose = FALSE))
cat(sprintf("%d workers: %.1f s\n", n, t[["elapsed"]]))
}
The sweeps above were produced by bench_tm_calculate.R,
driven by bench_tm_calculate_local.sh on the laptop and by
bench_tm_calculate.lsf on the compute node. All three are
in system.file("scripts", package = "TmCalculator"), and
their summaries are the two CSV files read here,
bench_hg38_laptop.csv and
bench_hg38_cluster.csv in extdata, so
re-running one on your own machine produces a file this vignette can
read.
bench_tm_calculate.R in the same directory times
tm_calculate() itself over a range of worker counts, with
bench_tm_calculate.lsf to submit it, and is the one to
reach for when sizing a new machine.
One function, two ways
of using it
tm_calculate() starts either from a source it can open,
as above, or from sequences you already hold:
tm_calculate(c("ACGTGCTAGCTAGCTAGC", "GGCCATATATGCGC"), method = "tm_nn", Na = 50)
Given sequences and nothing else it does exactly what it has always
done, one Tm per sequence, by the shortest path through the function.
Add regions, window or BPPARAM
and the profiling machinery engages.
What it will not do is divide the sequences of one region among
workers. With the compiled core the per-window loop is a minority of a
call’s cost; window construction, sequence retrieval and result assembly
run once, and sending the sequences to workers costs more than the loop
it divides. Measured on chromosome 1, splitting one call across five
workers was never faster than not splitting it. Parallelism therefore
divides by region, and sequences handed in directly are staged to a
temporary file so that the workers read them rather than receive
them.