Skip to contents

BamScale

Bioconductor-friendly multithreaded BAM processing

Bioconductor Bioc build Bioc downloads R-CMD-check pkgdown License: MIT

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

  • 2.3–3.2× faster single-file reads than scanBam / readGAlignments
  • Drop-in: same objects, same ScanBamParam filtering, same BiocParallel model
  • 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

Benchmarks

Intel Xeon Gold 6252 (96 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 BamScale configuration vs the single-threaded standard reader:

Read pattern Comparator Standard BamScale Threads Speedup
Core alignment fields Rsamtools::scanBam 15.6 s 6.8 s 48 2.3×
GAlignments object GenomicAlignments::readGAlignments 10.9 s 3.4 s 48 3.2×
Sequence + base quality Rsamtools::scanBam 22.3 s 8.1 s 24 2.8×

End-to-end workflows — BamScale 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 93% 3.7×
Coverage → RleList 76% 2.5×
Coverage → bigWig 21% 1.2×

Every result is measured against output verified equal to the standard tool: the ATAC fragment-size table is byte-identical to ATACseqQC::fragSizeDist across 49.8M reads, and the coverage RleList is identical() to GenomicAlignments::coverage() at every thread count.

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

BamScale is on Bioconductor (currently in the devel branch):

if (!requireNamespace("BiocManager", quietly = TRUE)) install.packages("BiocManager")
BiocManager::install(version = "devel")
BiocManager::install("BamScale")

Or install the development version from GitHub:

remotes::install_github("cparsania/BamScale")

Quick start 

library(BamScale)

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)

Output modes & compatibility

BamScale 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

BamScale 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, BamScale 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 BamScale 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.

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 and ATAC QC).

Citation 

citation("BamScale")

If BamScale 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 BamScale_*.tar.gz
Rscript -e "BiocCheck::BiocCheck('.')"

Community & support

License

MIT — see LICENSE.