Bioconductor-friendly multithreaded BAM processing
BamScaleR is a multithreaded BAM reader for R, built on the ompBAM OpenMP engine. It returns the exact Bioconductor objects you already use — verified byte-identical to Rsamtools and GenomicAlignments — but decodes each file across many cores, taking the read step off the critical path in alignment-centric workflows.
-orange){kind=icon}
scanBam / readGAlignments, and never slower — it matches or beats them even on a single threadScanBamParam filtering, same BiocParallel modelfragment_sizes(), mapq_dist(), bam_coverage(), and single-pass bam_coverage_bigwig() fold the computation inside the reader — byte-identical results without materialising per-read objectsthreads within a file, BiocParallel across filesA note on the name. This package was formerly called BamScale and was renamed to BamScaleR to avoid confusion with BAMscale (Pongor et al., Epigenetics & Chromatin 2020), an unrelated NCBI command-line tool that generates normalised coverage tracks and quantifies peaks. The R suffix marks this as the R/Bioconductor implementation; the two projects are not related.
Measured on a deeply sequenced public ENCODE ATAC-seq BAM (ENCFF191LNL, 151,057,932 alignments, 6.8 GB). Intel Xeon Gold 6252 (96 logical cores), warm page cache, median of 5 iterations. Full methodology, figures, and fair-comparison details are in the benchmark article.
Single-file read throughput — best BamScaleR configuration vs the single-threaded standard reader:
| Read pattern | Comparator | Standard | BamScaleR | Threads | Speedup |
|---|---|---|---|---|---|
| Core alignment fields | Rsamtools::scanBam |
263.3 s | 136.0 s | 36 | 1.94× |
GAlignments object |
GenomicAlignments::readGAlignments |
290.1 s | 115.9 s | 72 | 2.50× |
| Sequence + base quality | Rsamtools::scanBam |
409.8 s | 179.3 s | 24 | 2.29× |
End-to-end workflows — BamScaleR swapped in for the read step only, every other step identical on both sides. The gain tracks how read-bound the workflow is (Amdahl's law):
| Workflow | Read fraction | End-to-end speedup |
|---|---|---|
| ATAC fragment-size QC | 94% | 4.39× |
Coverage → RleList
|
72% | 2.49× |
| Coverage → bigWig | 41% | 1.52× |
Every result is measured against output verified equal to the standard tool: the ATAC fragment-size table is byte-identical to ATACseqQC::fragSizeDist across all 151,057,932 alignments, and the coverage RleList is identical() to GenomicAlignments::coverage() at every thread count.
In-reader aggregation — the reduction is folded into the multithreaded decoder, so no per-read R objects are built at all:
| Endpoint | Standard Bioconductor route | Standard | BamScaleR | Speedup |
|---|---|---|---|---|
| Fragment sizes | table(abs(scanBam(isize))) |
129.8 s | 20.9 s | 6.21× |
| MAPQ distribution | table(scanBam(mapq)) |
130.7 s | 23.9 s | 5.47× |
Coverage RleList
|
coverage(readGAlignments()) |
234.2 s | 33.2 s | 7.06× |
| BAM → bigWig | export.bw(coverage(readGAlignments())) |
410.2 s | 69.0 s | 5.95× |
bam_coverage_bigwig() also outperforms megadepth, a dedicated C++ coverage engine, by 3.45× at matched core counts (69.0 s vs 238.0 s at 72 threads) while leaving the track on disk and control in the R session.
Memory scales with the alignments present, not with the genome the BAM header declares: bam_coverage() on a 1.4 MB BAM carrying a whole-genome header needs 0.55 GB, not the 12.6 GB a header-sized allocation would require.
BamScaleR is not a raw-decode replacement for command-line
samtools. Its job is fast, object-faithful decoding inside R — putting the cores a single-threaded reader leaves idle to work, and handing back the precise Bioconductor objects your analysis depends on.
Requires R with a C++17 toolchain and an OpenMP-capable compiler; ompBAM and the other dependencies are pulled in automatically.
Install from GitHub:
remotes::install_github("cparsania/BamScaleR")
BamScaleR is under Bioconductor review (BiocContributions#226). Once it is accepted it will install with:
if (!requireNamespace("BiocManager", quietly = TRUE)) install.packages("BiocManager")
BiocManager::install(version = "devel")
BiocManager::install("BamScaleR")
library(BamScaleR)
bam <- ompBAM::example_BAM("Unsorted")
# Core alignment fields (scanBam-style extraction)
x <- bam_read(
bam,
what = c("qname", "flag", "rname", "pos", "mapq", "cigar"),
threads = 4
)
# A GAlignments object (drop-in for readGAlignments)
ga <- bam_read(
bam,
what = c("rname", "pos", "cigar", "strand"),
as = "GAlignments",
threads = 4
)
# Sequence + base quality, in standard-compatible objects
sq <- bam_read(
bam,
what = c("qname", "seq", "qual"),
seqqual_mode = "compatible",
threads = 4
)
# Fast chromosome-level counts
cnt <- bam_count(bam, threads = 4)
# In-reader aggregation: summaries computed inside the C++ reader,
# byte-identical to the scanBam/coverage equivalents
fs <- fragment_sizes(bam, threads = 4) # paired-end fragment-size distribution
mq <- mapq_dist(bam, threads = 4) # mapping-quality distribution
cov <- bam_coverage(bam, threads = 4) # per-base coverage RleList
bam_coverage_bigwig(bam, "coverage.bw", threads = 4) # single-pass BAM -> bigWig
BamScaleR fits where Rsamtools and GenomicAlignments already fit, adding a within-file threading axis on top:
what = c("qname", "flag", "rname", "pos", "mapq", "cigar", ...) — the metadata pattern behind most QC, filtering, and fragment-level summaries.as = "GAlignments" or as = "GAlignmentPairs".scanBam-shaped lists via as = "scanBam", or a data.frame / S4Vectors::DataFrame.ScanBamParam-style filtering: mapqFilter, flag, which, what, tag.BiocParallel (BPPARAM), exactly as in a standard workflow.BamScaleR parallelizes on two axes — across files via BPPARAM workers, and within each file via OpenMP threads. Approximate effective concurrency:
min(length(files), bpnworkers(BPPARAM)) * threads
With auto_threads = TRUE, BamScaleR keeps per-file thread counts high when possible — reducing the number of concurrently active file workers before it shrinks per-file threads. Within-file threading is the key difference from the baseline model, and the reason BamScaleR removes the read bottleneck when there are fewer files than cores.
compatible vs compact
seqqual_mode = "compatible" (default for equivalence) returns ordinary DNAStringSet / PhredQuality-style output, byte-identical to scanBam.
seqqual_mode = "compact" is a lower-level, deferred-decoding representation for when extraction throughput matters more than immediate string materialization:
seq — a list-column of raw vectors holding BAM-native packed bytes (two bases per byte)qual — a list-column of raw vectors holding per-base Phred bytes (255 = missing)qwidth is required to decode compact seq back to base lettersDecode it explicitly when needed:
sq_compact <- bam_read(
bam,
what = c("qname", "qwidth", "seq", "qual"),
seqqual_mode = "compact",
threads = 4
)
decoded <- decode_seqqual_compact(sq_compact)
param$which is implemented as sequential filtering, not indexed random-access jumps.seqqual_mode = "compact" is optimization-oriented and not intended for strict cross-package output equivalence.GAlignments / GAlignmentPairs outputs exclude unmapped records by design.samtools reach several-fold higher raw throughput, because they never construct the Bioconductor objects an analysis needs.Benchmark drivers and reporting assets ship under inst/benchmarks/ — see its README. Reported results come from run_server_benchmark.R (read patterns) and run_workflow_benchmark.R (end-to-end coverage, ATAC QC and the in-reader aggregation endpoints). The numbers above are from the run of 2026-10-08 on the 25-file ENCODE corpus, which is reconstructible from the MD5-verified manifest via download_atac_data.R.
citation("BamScaleR")
If BamScaleR contributes to a performance claim, please also cite ompBAM.
Pull requests are welcome — please include a short motivation, tests for behavior changes, and benchmark evidence for performance claims. Before opening a PR:
R CMD build .
R CMD check --as-cran BamScaleR_*.tar.gz
Rscript -e "BiocCheck::BiocCheck('.')"
MIT — see LICENSE.