Skip to contents

BamScaleR

Bioconductor-friendly multithreaded BAM processing

GitHub version Bioc devel version Bioc release R-CMD-check Bioc build pkgdown Bioconductor availability Bioc downloads License: MIT

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.

  • 1.9–2.5× faster single-file reads than scanBam / readGAlignments, and never slower — it matches or beats them even on a single thread
  • Drop-in: same objects, same ScanBamParam filtering, same BiocParallel model
  • In-reader aggregation: fragment_sizes(), mapq_dist(), bam_coverage(), and single-pass bam_coverage_bigwig() fold the computation inside the reader — byte-identical results without materialising per-read objects
  • Two parallel axes: OpenMP threads within a file, BiocParallel across files
  • Verified correct: output is byte-identical to the standard tools at every thread count

A 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.

Benchmarks

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.

Installation

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")

Quick start 

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

Output modes & compatibility

BamScaleR fits where Rsamtools and GenomicAlignments already fit, adding a within-file threading axis on top:

  • Field extraction via what = c("qname", "flag", "rname", "pos", "mapq", "cigar", ...) — the metadata pattern behind most QC, filtering, and fragment-level summaries.
  • Alignment objects via 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.
  • File-level parallelism via BiocParallel (BPPARAM), exactly as in a standard workflow.

Parallelism model

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.

Sequence/quality: 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 letters

Decode 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)

Current limitations

  • 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.
  • Intra-file thread scaling is strongly sublinear and saturates well below the full core count; the benefit is putting otherwise-idle cores to work on a single large file, not near-linear speedup.
  • Object-free command-line decoders such as samtools reach several-fold higher raw throughput, because they never construct the Bioconductor objects an analysis needs.

Reproducing the benchmarks

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 

citation("BamScaleR")

If BamScaleR contributes to a performance claim, please also cite ompBAM.

Contributing

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('.')"

Community & support

License

MIT — see LICENSE.