Benchmarking BamScale: reader throughput and end-to-end workflow impact
30 July 2026
Contents
0.1 Overview
BamScale is a multithreaded BAM reader, built on the ompBAM OpenMP engine, that
returns the same native Bioconductor objects as the standard readers and is a
drop-in replacement for Rsamtools::scanBam and
GenomicAlignments::readGAlignments.
Reading a BAM is only one step of an analysis, so a faster reader speeds a workflow end-to-end only in proportion to the time that workflow spends reading (Amdahl’s law). This vignette summarises benchmarks that quantify both the raw read speedup and its end-to-end effect on two standard workflows — coverage/bigWig track generation and ATAC-seq fragment-size QC — with BamScale output verified byte-identical to the standard tools.
The numbers below are a concise summary of a representative benchmark run
(Intel Xeon Gold 6252, 96 cores; warm page cache; median of 5 iterations). The
benchmark harness that produces them is in inst/benchmarks/
(run_server_benchmark.R, run_workflow_benchmark.R); see the last section to
reproduce the full results.
0.2 Read throughput on a single large BAM
On a single large BAM, BamScale reads 2.3–3.2x faster than the standard reader across three representative access patterns, because it spreads the BGZF decode across cores:
| Workload | Reads | Standard (s) | BamScale (s) | Threads | Speedup |
|---|---|---|---|---|---|
| Core fields | qname/flag/rname/pos/mapq/cigar | 15.6 | 6.8 | 48 | 2.29x |
| GAlignments | -> GAlignments object | 10.9 | 3.4 | 48 | 3.22x |
| Sequence + quality | seq + base quality | 22.3 | 8.1 | 24 | 2.76x |
Two honest points. At one thread BamScale is at or near the single-threaded
readers (the win comes from threading, not a faster single-core path), and the
scaling is strongly sublinear – roughly 2–3x on 48 threads, saturating by
about 24 threads. The value is not near-linear scaling but the ability to use cores
that a single-threaded reader cannot use at all – as on one large BAM in an
interactive session. Object-free command-line decoders such as samtools reach
several-fold higher raw throughput, but return no Bioconductor object; BamScale’s
role is fast, object-faithful decoding inside R.
0.3 End-to-end impact tracks the read fraction
Embedding BamScale as the read step of a complete workflow, the end-to-end speedup follows how read-dominated that workflow is (Amdahl’s law):
| Endpoint | Layer | Read fraction | Speedup |
|---|---|---|---|
| ATAC fragment-size QC | workflow | 93% | 3.66x |
| GAlignments (read) | read | 100% | 3.22x |
| Sequence + quality (read) | read | 100% | 2.76x |
| Coverage -> RleList | workflow | 76% | 2.50x |
| Core fields (read) | read | 100% | 2.29x |
| Coverage -> bigWig | workflow | 21% | 1.23x |
- ATAC fragment-size QC is nearly pure read, so BamScale’s read speedup carries through almost intact: 3.66x end-to-end.
- Coverage -> bigWig is dominated by the shared, single-threaded bigWig writer
(identical work on both arms), so the end-to-end gain is bounded at 1.23x even
though the read phase itself is much faster. Stopping at the in-memory coverage
RleList– what many analyses consume next – recovers a 2.50x gain.
0.4 A note on cores and fairness
A multithreaded reader can occupy cores a single-threaded one cannot. At matched core counts (both arms given the same number of cores) the two readers are approximately at parity for multi-file processing (coverage ~1.10x, ATAC ~0.96x): the single-threaded reader already saturates a machine by parallelising across files, so BamScale’s multi-file advantage comes specifically from also threading within a file – which pays off when there are fewer files than cores.
0.5 Output is identical to the standard tools
Every speedup above is against output verified equal to the standard tool. The
coverage RleList is identical() to GenomicAlignments::coverage() on
readGAlignments output at every thread count, and the ATAC fragment-size table is
identical() to Rsamtools::scanBam and byte-identical (maximum absolute
difference 0) to ATACseqQC::fragSizeDist over 49.8 million reads. This is what
makes BamScale a genuine drop-in rather than an approximation.
0.6 Try it
library(BamScale)
bam <- ompBAM::example_BAM("Unsorted")
## Alignment fields as a GAlignments object (drop-in for readGAlignments)
ga <- bam_read(bam, what = c("rname", "pos", "cigar", "strand"),
as = "GAlignments", threads = 2)
ga## GAlignments object with 10000 alignments and 0 metadata columns:
## seqnames strand cigar qwidth start end
## <Rle> <Rle> <character> <integer> <integer> <integer>
## [1] 19 + 1S88M6798N61M 150 572614 579560
## [2] 19 - 1S148M1S 150 579499 579646
## [3] 6 + 60S86M4S 150 44252112 44252197
## [4] 6 - 1S146M3S 150 44252124 44252269
## [5] 12 + 1S149M 150 46185884 46186032
## ... ... ... ... ... ... ...
## [9996] 2 - 16S122M12S 150 186680457 186680578
## [9997] 15 + 4S143M3S 150 90501177 90501319
## [9998] 15 - 26M2I56M66S 150 90501301 90501382
## [9999] 4 + 3M336N95M752N51M1S 150 121801488 121802724
## [10000] 4 - 1S117M121N32M 150 121802677 121802946
## width njunc
## <integer> <integer>
## [1] 6947 1
## [2] 148 0
## [3] 86 0
## [4] 146 0
## [5] 149 0
## ... ... ...
## [9996] 122 0
## [9997] 143 0
## [9998] 82 0
## [9999] 1237 2
## [10000] 270 1
## -------
## seqinfo: 25 sequences from an unspecified genome## ...or core fields as a data.frame; `threads` controls within-file parallelism
df <- bam_read(bam, what = c("qname", "flag", "mapq"),
as = "data.frame", threads = 2)
head(df)## qname flag mapq
## 1 ST-E00600:137:H77Y3CCXY:1:1101:6837:1309 163 255
## 2 ST-E00600:137:H77Y3CCXY:1:1101:6837:1309 83 255
## 3 ST-E00600:137:H77Y3CCXY:1:1101:10450:1309 163 255
## 4 ST-E00600:137:H77Y3CCXY:1:1101:10450:1309 83 255
## 5 ST-E00600:137:H77Y3CCXY:1:1101:15077:1309 99 255
## 6 ST-E00600:137:H77Y3CCXY:1:1101:15077:1309 147 2550.7 Reproducing the full benchmarks
The complete benchmark harness ships under inst/benchmarks/:
dir(system.file("benchmarks", package = "BamScale"))
# run_server_benchmark.R : read-pattern micro-benchmarks (step1 / GAlignments / seq+qual)
# run_workflow_benchmark.R : end-to-end coverage and ATAC fragment-size QC workflows
# download_atac_data.R : fetch + index the ENCODE GM12878 ATAC BAMs used hereEach writes machine-readable result tables (summary.csv) plus host and
correctness metadata, from which the summary figures above are derived.
0.8 Session information
## R version 4.6.1 (2026-06-24)
## Platform: x86_64-pc-linux-gnu
## Running under: Ubuntu 24.04.4 LTS
##
## Matrix products: default
## BLAS: /home/biocbuild/bbs-3.24-bioc/R/lib/libRblas.so
## LAPACK: /usr/lib/x86_64-linux-gnu/lapack/liblapack.so.3.12.0 LAPACK version 3.12.0
##
## locale:
## [1] LC_CTYPE=en_US.UTF-8 LC_NUMERIC=C
## [3] LC_TIME=en_GB LC_COLLATE=C
## [5] LC_MONETARY=en_US.UTF-8 LC_MESSAGES=en_US.UTF-8
## [7] LC_PAPER=en_US.UTF-8 LC_NAME=C
## [9] LC_ADDRESS=C LC_TELEPHONE=C
## [11] LC_MEASUREMENT=en_US.UTF-8 LC_IDENTIFICATION=C
##
## time zone: America/New_York
## tzcode source: system (glibc)
##
## attached base packages:
## [1] stats graphics grDevices utils datasets methods base
##
## other attached packages:
## [1] ggplot2_4.0.3 BamScale_0.99.13 BiocStyle_2.41.0
##
## loaded via a namespace (and not attached):
## [1] SummarizedExperiment_1.43.0 gtable_0.3.6
## [3] xfun_0.60 bslib_0.11.0
## [5] Biobase_2.73.2 lattice_0.22-9
## [7] vctrs_0.7.3 tools_4.6.1
## [9] bitops_1.1-0 generics_0.1.4
## [11] stats4_4.6.1 parallel_4.6.1
## [13] tibble_3.3.1 pkgconfig_2.0.3
## [15] Matrix_1.7-6 RColorBrewer_1.1-3
## [17] S7_0.2.2 S4Vectors_0.51.6
## [19] cigarillo_1.3.1 lifecycle_1.0.5
## [21] compiler_4.6.1 farver_2.1.2
## [23] Rsamtools_2.29.0 Biostrings_2.81.6
## [25] tinytex_0.60 Seqinfo_1.3.0
## [27] codetools_0.2-20 GenomeInfoDb_1.49.1
## [29] htmltools_0.5.9 sass_0.4.10
## [31] yaml_2.3.12 pillar_1.11.1
## [33] crayon_1.5.3 jquerylib_0.1.4
## [35] BiocParallel_1.47.0 DelayedArray_0.39.3
## [37] cachem_1.1.0 magick_2.9.1
## [39] abind_1.4-8 ompBAM_1.17.0
## [41] tidyselect_1.2.1 digest_0.6.39
## [43] dplyr_1.2.1 bookdown_0.47
## [45] labeling_0.4.3 fastmap_1.2.0
## [47] grid_4.6.1 cli_3.6.6
## [49] SparseArray_1.13.2 magrittr_2.0.5
## [51] S4Arrays_1.13.0 dichromat_2.0-1
## [53] withr_3.0.3 UCSC.utils_1.9.0
## [55] scales_1.4.0 rmarkdown_2.31
## [57] XVector_0.53.0 httr_1.4.8
## [59] matrixStats_1.5.0 otel_0.2.0
## [61] evaluate_1.0.5 knitr_1.51
## [63] GenomicRanges_1.65.1 IRanges_2.47.2
## [65] rlang_1.3.0 Rcpp_1.1.2
## [67] glue_1.8.1 BiocManager_1.30.27
## [69] BiocGenerics_0.59.10 jsonlite_2.0.0
## [71] R6_2.6.1 MatrixGenerics_1.25.0
## [73] GenomicAlignments_1.49.1