When not to reach for advanced libraries
Adapted from a tutorial I gave at a lab group meeting. The full workshop has the runnable code.
The previous post ended on a recommendation: store your data as Parquet and use an engine that can query it without loading it. On a grouped aggregation over 46 million rows, a DuckDB pipeline finished in half a second using 0.14 GB where data.table needed 9.9 seconds and 14 GB, and almost all of data.table’s disadvantage was the loading rather than the computing.
Here we explore four more standard bioinformatic types of workflows, and see how they perform with these different tools.
Every table below splits the time the same way as the previous post. Load is getting the inputs into the form the method computes on, which for the engines means little more than opening a handle. Compute is the actual work, and for the engines it includes reading the file, because that is inseparable from the query by design. Total additionally includes R starting up and attaching packages, which is a few tenths of a second in every row.
Same conditions throughout: 46 million somatic mutation calls, 2,778 samples, each method in its own fresh R process, with memory sampled from outside R. The script and its raw output are alongside this post.
Three tables feed the workflows. Variant calls, one row per mutation:
| samplename | chromosome | position | cluster_id |
|---|---|---|---|
| 65b0b887db6f… | 17 | 74655400 | 0 |
| 65b0b887db6f… | 5 | 114200192 | 0 |
| 65b0b887db6f… | 1 | 9923399 | 0 |
Gene models, 20,000 rows:
| gene_id | symbol | chromosome | gene_start | gene_end |
|---|---|---|---|---|
| ENSG00000000001 | GENE1 | 1 | 302322 | 420794 |
| ENSG00000000023 | GENE23 | 1 | 532631 | 652894 |
And a catalogue of known positions, 36.3 million rows:
| chromosome | pos | rsid | allele_freq |
|---|---|---|---|
| 1 | 26 | rs18955602 | 0.9678 |
| 1 | 87 | rs30707517 | 0.5598 |
Assigning variants to genes: the range join
You have a variant call set and want to know which gene each mutation falls in, so you can count mutations per gene and look for genes hit more often than chance would predict. That means labelling each of 46 million positions with the gene whose start and end bracket it. This is an interval overlap rather than an equality match, and it is one of the most common operations in genomics.
| Method | Load | Compute | Total | Peak RSS (GB) |
|---|---|---|---|---|
| data.table (binned key) | 1.33 | 1.34 | 3.12 | 2.25 |
data.table (foverlaps) |
1.38 | 6.51 | 8.36 | 2.23 |
| DuckDB (range join) | 0.05 | 8.27 | 8.67 | 0.23 |
GenomicRanges (findOverlaps) |
4.12 | 6.13 | 11.48 | 3.34 |
| Polars (join + filter) | 0.02 | 121.01 | 121.33 | 1.44 |
Those are five genuinely different strategies, which the names do not make obvious:
- DuckDB expresses the overlap as SQL
BETWEEN. A column store joins by scanning and hashing, and neither of those answers an inequality, so it evaluates the predicate against candidate pairs. foverlapsbuilds an interval tree over the 20,000 genes and probes it once per variant.findOverlapsis the Bioconductor equivalent, also an interval index.- The binned key cuts the genome into 10 kb windows, copies each gene into every window it spans, joins on chromosome and window as an ordinary equality join, then filters the survivors on the real coordinates.
- Polars has no interval structure available here, so it equi-joins on chromosome and filters afterwards.
Both interval-tree methods beat DuckDB on the compute step, 6.51 and 6.13 seconds against 8.27, so this is not an artifact of loading. A column store is fast because it scans and hashes large volumes quickly, and a range join does neither. Range queries are common in genomics, so this case comes up often.
GenomicRanges is the interesting row. It has the fastest compute of any method that is not binned, and still the worst total of the three, because it spends 4.12 seconds building GRanges objects first. The S4 machinery that makes it compose with the rest of Bioconductor is not free, but if you are already working within this environment then this is a cost you already have paid.
The memory column still favors DuckDB, which does the job in 0.23 GB against 2.23 to 3.34 for the in-memory approaches, because it never materializes the variants in R. On time these are close. On memory they are quite different, which favors the situation where you are trying to run this analysis in a limited environment.
Two other things in that table are worth noting.
Avoid bad approaches. Without a specialized function, one approach would be to join on chromosome and then filter on position. Joining 46 million variants to 20,000 genes on chromosome alone produces 41.9 billion intermediate rows, roughly 1.4 TB at 36 bytes each, before the filter removes almost all of them. That is the Polars row. Streaming is the only reason it finishes: it holds its peak to 1.44 GB and takes 121 seconds, about 15× the foverlaps total. Anything that materializes those pairs runs out of memory.
Reformulating the problem gains efficiency. The binned row turns the range join into an equality join and filters afterwards. Its compute step is 1.34 seconds against DuckDB’s 8.27, although DuckDB still wins out on memory.
A result with millions of rows
The workflow above returned one row per gene. Here, the output is similar in size to the input.
Example analysis: Mutation density along the genome. For each sample, count how many mutations fall in each 1 Mb window. That profile is what you plot as a mutation rate along a chromosome, or scan for hypermutated regions. The operation is a group-by on sample, chromosome, and position divided by a million:
| samplename | chromosome | bin | n |
|---|---|---|---|
| 000b42459531… | 1 | 2 | 1 |
| 000b42459531… | 1 | 5 | 2 |
| 000b42459531… | 1 | 6 | 1 |
| 000b42459531… | 1 | 8 | 1 |
46,115,588 rows go in and 8,870,305 come out. Both are large.
| Method | Load | Compute | Total | Peak RSS (GB) |
|---|---|---|---|---|
| DuckDB | 0.03 | 0.83 | 1.34 | 1.13 |
| Polars (streaming) | 0.00 | 1.93 | 2.34 | 2.67 |
| data.table | 3.93 | 1.14 | 5.59 | 4.81 |
This is where splitting the compute loads by phase highlights tradeoffs. On total time DuckDB wins by 4.2×. On compute the gap is 0.83 against 1.14, about 37%, and data.table is faster than Polars, while using the most memory.
Nearly all of DuckDB’s advantage is that data.table spent 3.93 seconds loading a file DuckDB never loaded. The compute gap is much smaller than in the previous post because a large share of the engine’s time now goes to serializing 8.9 million result rows into R’s memory rather than to aggregating.
The engine’s advantage is largest when the result is small, and it shrinks as the result grows, because the handover into R costs more. If a large intermediate is only going to feed another query, leave it in the engine or write it straight to Parquet with sink_parquet() or COPY ... TO.
Joining two large tables
You want to know which of your variants are already catalogued in something like dbSNP and which look novel, counted per sample. That means flagging which of the 46 million variants appear in a 36-million-row table of known positions, matched on chromosome and position together.
| Method | Load | Compute | Total | Peak RSS (GB) |
|---|---|---|---|---|
| DuckDB | 0.05 | 0.73 | 1.16 | 2.96 |
| Polars (streaming) | 0.00 | 2.30 | 2.61 | 4.13 |
data.table (merge) |
27.38 | 16.01 | 43.85 | 6.36 |
| data.table (packed key) | 35.69 | 9.33 | 45.48 | 6.37 |
base R (merge) |
3.26 | 316.35 | 320.08 | 23.06 |
Joining two large tables on a compound key is the case where an in-memory approach becomes impractical rather than just slow. Both sides have to be materialized first, which is why data.table spends 27 seconds loading: it reads 46 million variants, reads 36 million catalogue rows, and sorts the second to set a key. The join then builds a hash table, and the peak sits well above the sum of the inputs.
The two merge rows put the same cost in different columns. data.table sorts up front in setkey(), so the expense lands in the load column. Base R reads the same two files in 3.26 seconds and then spends 316 seconds inside merge(), doing a sort-based join with no index.
That last row peaked at 23.06 GB for a join whose inputs total under half a gigabyte on disk. On a 16 GB machine it would simply fail.
DuckDB and Polars stream both sides from Parquet, never hold either in full, and spill to disk if the hash table outgrows the budget.
The composite key
Joining on two columns invites a shortcut: combine them into one key and join on that instead.
That is the data.table (packed key) row, and it returns the same answer as every other method in the table.
It is 1.7× faster on the join itself, 9.33 seconds against 16.01, because hashing one numeric column is cheaper than hashing a pair. The whole gain then goes straight back out in the load column: building the key costs an extra pass over both tables and two new vectors of 46 and 36 million elements, which pushes loading from 27.38 to 35.69 and leaves the total slightly worse than the plain two-column join.
At this scale it is not worth doing. It pays off when you reuse the key across several joins, since the construction cost is paid once and the 1.7× applies every time.
A full analysis end to end
You want to test whether tumors with a high subclonal mutation fraction differ in the expression of a few genes you care about. That is an analysis rather than a single operation: split samples by subclonal mutation fraction (sMF), then compare expression of five key genes between the high and low groups. Counts per million require library sizes, which require summing all 55.5 million count rows, even though the question at the end concerns five genes.
| Method | Load | Compute | Total | Peak RSS (GB) |
|---|---|---|---|---|
| DuckDB (one query) | 0.03 | 0.22 | 0.63 | 0.21 |
| R in memory (naive) | 3.88 | 7.57 | 12.02 | 6.70 |
DuckDB is 19× faster on total time and uses 32× less memory for an identical answer. It wins on both phases here: 3.88 seconds of loading it never does, and 0.22 against 7.57 on the compute.
The naive pipeline reads all 55.5 million count rows into R because it needs the library sizes, joins them, then throws away everything except five genes. Its peak memory is set by the largest object it held along the way, which has nothing to do with the size of the answer.
A large computation feeding a narrow question is what a query engine handles best, and it is a common pattern in genomics.
Key takeaways
These workflows do not point at one rule. They separate into cases where the engines help a great deal, cases where they help a little, and cases where the choice does not matter much.
Engines help most when the operation is filtering, joining or aggregating, the input is larger than you want to hold in memory, and the answer is much smaller than the input. The end-to-end analysis is the clearest example, 19× faster and 32× lighter. The 46M by 36M join is the other, where the in-memory versions need 6 to 23 GB and the engines need 3 to 4.
Memory is a more reliable reason to switch than speed. Across these tables the engines are sometimes faster and sometimes slower, but they are consistently lighter, often by an order of magnitude. A job that is 30% slower still finishes, while a job that needs 23 GB on a 16 GB machine does not.
The advantage narrows as the result grows. Returning 8.9 million rows cut the compute gap to 37% while the total gap stayed above 4×, because by then the difference was loading rather than computing.
Some work has no relational form at all. Matrix operations are a good example.These are the operations that R is built to perform, and SQL is not a good fit.
If your data is already in memory, most of this does not apply. The load column is what separates the engines from the in-memory tools in nearly every table here. Working interactively on a data frame you read in an hour ago, that cost is already paid and the remaining differences are much smaller.