--- title: "Genome-Scale Tm Profiling on hg38: Inputs, Parallel Execution and Performance" author: "Junhui Li, Lihua Julie Zhu" date: "`r Sys.Date()`" output: rmarkdown::html_vignette: toc: true toc_depth: 3 number_sections: true vignette: > %\VignetteIndexEntry{Genome-Scale Tm Profiling on hg38: Inputs, Parallel Execution and Performance} %\VignetteEngine{knitr::rmarkdown} %\VignetteEncoding{UTF-8} --- ```{r setup, include=FALSE} knitr::opts_chunk$set( echo = TRUE, eval = FALSE, # the hg38 chunks need BSgenome.Hsapiens.UCSC.hg38 and minutes message = FALSE, warning = FALSE, fig.retina = 1, dpi = 72 ) ``` # 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: 1. **What can I give it?** A BSgenome package, a FASTA file, or sequences already in R, and any set of regions within them. 2. **How does it divide the work?** Into tasks, one region or one segment each, dispatched to workers. 3. **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. ```{r install, eval=FALSE} 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: ```{r quick-start} 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: ```{r quick-start-parallel} 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. ```{r sources} 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. ```{r regions} ## 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. ```{r unit} 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. ```{r sweep-data, eval=TRUE} # 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), ] ``` ```{r sweep-table, eval=TRUE} 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.")) ``` ```{r sweep-figure, eval=TRUE, fig.width=8, fig.height=3.6, fig.cap="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."} 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") } par(op) ``` 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: ```{r sweep-your-machine} 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. # What comes back A `GRanges` with `Tm` and `GC` metadata columns in genomic order, ready for `integrate_granges()`, `compare_groups()` and the plotting functions, and for anything else that takes genomic intervals. ```{r downstream} tm_annot <- integrate_granges(gr_tm = tm, gr_features = atac_peaks, strategy = "overlap", weight = "overlap") compare_groups(tm_annot, value_cols = "Tm", group_col = "class") ``` The sequence and complement columns are dropped by default, since they run to roughly 500 MB per large chromosome; pass `keep_sequence = TRUE` if you need them. # 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: ```{r tm-calculate} 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. # Session information ```{r session-info, eval=TRUE} sessionInfo() ```